I 


NASA 


Cs4 

CO 

CO 


I 



l/J 



LOAN COPY: RETURN TO 4 
AFWL (WLIL-2) 

kirtland afb, n mex 


A METHOD FOR ANALYZING 
THE AEROELASTIC STABILITY 
OF A HELICOPTER ROTOR 
IN FORWARD FLIGHT 

by Peter Crimi 

Prepared by 

ROCHESTER APPLIED SCIENCE ASSOCIATES, INC. 

Rochester, N. Y. 

for Langley Research Center 




NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


WASHINGTON, D. C. 


AUGUST 1969 


TECH LIBRARY KAFB, NM 



TECH LIBRARY KAFB, NM 


II 


□DbD505 


NASA. CR-1332 


A METHOD FOR ANALYZING THE AEROELASTIC STABILITY 
OF A HELICOPTER ROTOR IN FORWARD FLIGHT 


By Peter Crimi 


Distribution of this report is provided in the interest of 
information exchange. Responsibility for the contents 
resides in the author or organization that prepared it. 


Issued by Originator as RASA Report No. 68-10 


Prepared under Contract No. NAS 1-7411 by 
ROCHESTER APPLIED SCIENCE ASSOCIATES, INC. 

Rochester, N.Y. 

for Langley Research Center 
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 




CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION 2 

SYMBOLS 3 

DEVELOPMENT OF THE EQUATIONS OF MOTION OF A ROTOR IN 

FORWARD FLIGHT 6 

DEVELOPMENT OF THE BASIC METHOD 17 

OUTLINE OF BASIC COMPUTER PROGRAM 26 

APPLICATION OF THE METHOD 31 

CONCLUDING REMARKS 42 

APPENDIX A: Relationship Among the Solutions of 

Original and Differentiated Systems 4 3 

APPENDIX B: Listing of Basic Computer Program 4 5 

APPENDIX C: Method for Extraction of the Common Polynomial 

Factor From the Two Higher-Degree 

Polynomials 10 1 

REFERENCES 106 


iii 


A METHOD FOR ANALYZING THE AEROELASTIC 


STABILITY OF A HELICOPTER ROTOR IN 
FORWARD FLIGHT 

By Peter Crimi 

ROCHESTER APPLIED SCIENCE ASSOCIATES, INC. 

SUMMARY 


A method has been developed which provides the exact solution 
to the perturbation equations of motion of a rotor in forward 
flight. Effects of compressibility, stall and reversed flow are 
taken into account, and there is no restriction on the number of 
degrees of freedom which can be analyzed. 

The equations of motion of a rotor blade in forward flight 
were derived in terms of the normal modes of free vibration of the 
blade. Aerodynamic forces were expressed on the basis of quasi- 
steady flow, using strip theory. Perturbation equations were gen- 
erated by expanding all aerodynamic forces in the perturbation 
variables about their steady-state values. The equations are of 
the form of a coupled set of linear, second-order differential 
equations with periodically varying coefficients. 

A method was then developed for analyzing the stability of 
linear dynamic systems with periodically varying parameters. 
Stability is determined by calculation of all the characteristic 
values of the system. Thus, a solution provides rates of growth 
or decay of the motion following a disturbance, in addition to a 
determination of whether or not the system is stable. 

A digital computer program was prepared to implement the gen- 
eral method for the case of a rotor blade with one, two or three 
degrees of freedom. Stability calculations were performed for 
comparison with results obtained by direct time integration of the 
nonlinear equations of motion of a rigid blade with flapping and 
lead-lag hinges. Limited calculations were also made of the aero- 
elastic stability of a model rotor blade for which experimental 
flutter data was available. Comparisons of the results indicate 
qualitative agreement in both cases. 



INTRODUCTION 


The aeroelastic stability of helicopter rotors is of concern 
for two reasons. First, at the high advance ratios associated 
with compound and stowable-rotor operation, the hostile aerodynamic 
environment could lead to dangerous instabilities. Secondly, in 
conventional operation, instabilities can occur which, while gen- 
erally not of a catastrophic nature because of nonlinear effects, 
are nonetheless serious from a fatigue and control standpoint. 

The analysis of rotor flutter for hovering flight is a rela- 
tively straightforward problem, the formulations being essentially 
the same as those for classical flutter of conventional aircraft. 
Considerable research, both theoretical and experimental, has dealt 
with the flutter of hovering rotors (see, for example. References 1 
through 4) . Agreement between theory and experiment has generally 
been good. 

The problem for a rotor in forward flight is fundamentally 
different from that of the hovering case, due to the nature of the 
equations of motion. Because the relative flow imposed on a given 
blade section varies periodically as the blade rotates, the effec- 
tive dynamic pressure, and hence the constant of proportionality of 
the aerodynamic forces, also varies periodically. As a consequence, 
the equations of motion of the rotor have periodically varying co- 
efficients. Systems described by equations of this type, even though 
linear, display many unusual properties and in some respects resem- 
ble nonlinear systems (see Reference 5) . 

Linear differential equations with periodically varying coeffi- 
cients have been the concern of applied mathematicians for over a 
century. The differential equation of second order bearing his 
name was discussed by Mathieu in 1868 in reference to the problem 
of a vibrating elliptic membrane. The more general second-order 
equation derived by Hill for determining the motion of the lunar 
perigee was presented by him in 1877. In 1883, Floquet determined 
the form of the solution for any linear differential equation with 
periodic coefficients.* 

More recent analyses related to dynamic systems with periodic 
parameters are presented in References 7, 8 and 9. In Reference 

th 

7, the general N— - order problem is treated, while Reference 8 is 
concerned with spin-stabilized satellites and Reference 9 deals 

*A discussion of’ the” early history of MathieiT and related 
functions is given in Reference 5. Many papers of historical 
interest are also noted in Reference 6. 
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with mechanical instabilities of helicopter rotors. These analyses 
and similar ones treating related problems generally either rely on 
expansion in a small parameter to obtain a solution or are re- 
stricted to special forms of the differential equations. 

Since a general method for analyzing linear systems with peri- 
odic parameters has not been available, it has been necessary in the 
past, in investigating rotor flutter, to rely primarily on direct 
integration of the equations of motion using an analog or digital 
computer. The results of analog studies of rotor stability are re- 
ported in References 10 and 11. Digital computers were used to 
obtain results discussed in References 12 and 13. Computer solu- 
tions are useful particularly in that nonlinear effects can readily 
be incorporated to give an accurate indication of system dynamics 
for a specific rotor .and flight conditions. However, the results 
are often difficult to interpret and evaluate in terms of system 
stability, and, especially when a digital computer is used, their 
cost precludes conducting a thorough parametric study. 

It has been possible to gain some insight into rotor aero- 
elastic characteristics by treating a single degree of freedom 
analytically (References 14 and 15) . A single second-order differ- 
ential equation represents the system. This equation can be reduced 
to the form of Hill's equation (Reference 6) , for which there are 
techniques available to obtain solutions. The analyses have been 
confined to the flapping degree of freedom; it was found that this 
simple system is quite stable, but that flutter can occur at very 
high advance ratios (of the order of unity) . 

The study reported here was directed to obtaining the exact 
analytical solutions of the perturbation equations of motion of a 
rotor in forward flight. The values of system parameters or flight 
conditions were not restricted, nor were limitations placed on the 
number or types of degrees of freedom. Because the general problem 
was attacked, the solution obtained has application in areas other 
than rotor aeroelastic stability. The method presented can be used 
to determine the stability of any linear dynamic system with peri- 
odically varying parameters. Also, the computer program developed 
to implement the method has been segmented in such a way that it 
can be applied to the stability analysis of any linear dynamic 
system with three degrees of freedom. 


SYMBOLS 


a mn ' ^mn 


dimensionless, periodic coefficients in the basic 
set of differential equations 


C / 

mn 


mn 


dimensionless, periodic coefficients in the 
derived set of differential equations 
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blade chord, m 

blade section lift coefficient 
blade section drag coefficient 

blade section moment coefficient for moment about 
mid-chord, positive to increase incidence 

drag per unit span of blade, N/m 

aerodynamic force component per unit blade span 
in the y-direction, N/m 

aerodynamic force component per unit blade span 
in the z-direction, N/m 

coefficients in the finite sum of trigonometric 
functions related to A 

mass moment of inertia about the elastic axis 
per unit span; kg-m 

lift per unit span of blade, N/m 

horizontal distance forward of the X 0 axis of the 
elastic axis (x-axis) , m 

aerodynamic moment per unit span about mid-chord, 
positive to increase incidence, N 

aerodynamic moment per unit span about the elastic 
axis, positive to increase incidence, N 

blade mass per unit span, kg/m 

number of degrees of freedom of the elasto- 
mechanical system 

rotor radius, m 

magnitude of resultant fluid velocity relative 
to a blade section, m/s 

component of fluid velocity, relative to a blade 
section, normal to the Xo-Zq plane, m/s 

component of fluid velocity, relative to a blade 
section, tangent to the Xq-Zq plane, m/s 



magnitude of free-stream velocity (aircraft 
forward speed) , m/s 

displacement of blade section in the y-direction ,m 
displacement of blade section in the z-direction,m 
wake-induced inflow, m/s 

coordinates rotating with the blade, the Yo axis 
coincident with the shaft 

coordinates fixed with respect to a local blade 
section, the x-axis coinciding with the elastic 
axis and the z-axis parallel to the chord 

distance of the elastic axis ahead of mid-chord, m 

blade section angle of attack, rad 

shaft tilt with respect to normal to free stream, 
positive for aft tilt of the negative Yo-axis, rad 

infinite determinant 

distance forward of the elastic axis of blade 
section mass center, m 

displacement of the n— coupled mode of free 
vibration of the elasto-mechanical system 

flapwise bending slope 

real part of system characteristic value iw 
imaginary part of system characteristic value iw 


roots defining points at which A is singular 

mass density of free stream, kg/m 3 

angle between local blade chord and the Xo~Yo 
plane, rad 

torsional deflection about the elastic axis, posi- 
tive to increase incidence, rad 

chordwise bending slope 

rotor rotational speed, rad/s 



3 . 0 ) 


characteristic value of the basic or derived 

systems, iio = \ + ix 

JK X 

natural frequency of the k— coupled mode of 
free vibration of the rotating system, rad/s 


DEVELOPMENT OF THE EQUATIONS OF MOTION 
OF A ROTOR IN FORWARD FLIGHT 


The perturbation equations of motion for a rotor blade in 
forward flight are derived in this section. The set of coupled, 
linear differential equations with periodic coefficients which are 
obtained form the subject for analysis in the next section. 

The elasto-mechanical segment of the aeroelastic system, 
assumed to have N degrees of freedom, is conveniently represented 
in terms of the N normal modes of free vibration. Normal modes 
and frequencies for a rotating beam can be obtained either from a 
continuous representation (Reference 16) or from a lumped-mass 
model (Reference 17) . 

Aerodynamic forces are derived here in accordance with the 
usual strip-theory assumptions. Quasi-steady flow is assumed as 
well. A perturbation analysis is then performed, with expansion 
of experimentally or otherwise determined aerodynamic coefficients 
in a Taylor series, about the nominal flow conditions, in the quan- 
tities defining the perturbations in the blade motions. 

Consider a rotor blade, then, with angular speed fi subjected 
to a uniform free stream of magnitude and density p, as shown in 

Figure 1. The X 0 -axis is taken coincident with the pitch, or 

feathering, axis, when the blade is in its undeformed position, as 
shown in the figure. The coordinates (x,y,z) also rotate with the 
blade and are fixed with respect to a local blade section of the 
undeformed blade, with origin at the elastic axis and with the 
z-axis parallel to the chordline. The horizontal offset of the 
elastic axis with respect to the Xg-axis is denoted The angle 

which a section of the undeformed blade makes with the X 0 -Z 0 plane, 

combining built-in twist and pitch control settings, is denoted $. 


Blade motions are defined by the five variables v, w, <j> , 6 and 
i h The deflections of the elastic axis in the y and z directions 
are v and w, respectively, the torsional deflection is <j>, while 0 
and i p are bending slopes : 
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These five variables are expressed in terms of the N coupled modes 
of free vibration selected to represent the rotating mechanical 
system. Specifically, these variables are expressible in the form 


v(x,t) 


w(x,t) 


4> (x,t) 


e (x,t) 


<l> (x,t) 


l A v (n) <*> 5 (t) 
n=i 

! A^ n) (x)C n (t) 
n= i 

! »e n) (X)t n (t) 

n= i 


<D 


where (x) , (x) , etc. , denote the values of v, w, etc. , 

V ** 

obtained on excitation of the n— coupled mode. 

If Lagrange's equations are applied to the system, with the 
functions Si, £2/ •••/ S N as generalized coordinates, the system 

equations of motion are given by 


t w k — (t) , k - 1 ,2 , . . . ,N; (2) 


where is the natural frequency of the k— mode and is the 

generalized force applied to that mode. It is assumed that the 
mode shapes have been normalized so as to yield a generalized 
mass of unity. 
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The generalized forces are expressed by 

f R 

F k (t) = {A^, k * (x)Q v (x,t) + A^ k) (x)Q w (x,t) 

« 

0 

+ A^ k * (x)Q^(x,t) + A^ k ^ (x) (x,t) }dx 

(3) 

where R is rotor radius. Omitting the steady-state forcing terms, 
which are not of concern to a stability analysis, the force and 
moment distributions Q , Q , etc. are given by 


Q v (x,t) = - 2m£2e sin $ip + F v (x,t) 

Q (x,t) = - 2mQe cos $>iji + F (x,t) 
w w 

1 (4) 

Q (x , t) = 2I 0 fi sin l + M, (x,t) 

V <p 

Q^(x,t) = - 2l 0 n4i + 2mfJe [v sin 0 + w cos $] 

where F v and F w are aerodynamic forces in the y and z directions, 

respectively, is the aerodynamic moment about the elastic axis, 

e is the distance of the section mass center forward of the elastic 
axis, m is blade mass per unit span and I Q is mass moment of _ inertia 
about the elastic axis per unit span. The terms containing ip and $ 
represent gyroscopic coupling terms. They are placed on the right- 
hand side because their inclusion in the vibration analysis to 
obtain mode shapes and frequencies would disrupt the orthogonality 
of the modes. No quantity Q , similar to the quantity Q , appears 

in Eqs. (3) or (4) because the mass moment of inertia of a blade 
section about its chordline has been assumed to be negligible. 

The next step is to obtain expressions for the aerodynamic 
forces and to perform the appropriate linearizations. To accom- 
plish this define, first, an effective section incidence a. 


- 1 

a == $ + <j) + tan 
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where V and V. are components of fluid velocity relative to the 
n t 

blade, at midchord, normal and tangential to the rotor plane, 

respectively. If a denotes shaft angle (positive for aft tilt of 

s 

the D-vector) , Z is distance of the elastic axis forward of mid- 

cl 

chord and ftt is the blade azimuth relative to the downstream direc- 
tion, then these velocity components are given by 


V 


n 


[v + Z 4> - ft (Z - £ ) ( 0 cos $ + ^ sin $) ] cos $ 
a a z 

- w sin $ + Vj cos a g cos fit(e cos $ + i|) sin <f) 

+ V_ sin a - w. 
r si 

fix + V.c cos a sin fit + w cos $ 
f s 

+ [v + Z - fi(Z - l ) (0 cos $ + i|> sin <f) ] sin <5 
a a z 


} ( 5 ) 


where w^ is the wake-induced inflow, assumed constant over the 

rotor plane. The factor fi(e cos $ + ip sin $) appears in Eqs. ( 5 ) 
because the angular velocity directed along the Yq (shaft) axis, 
of magnitude fi, has a component which is tangent to the blade 
when the blade is inclined to the Xq - Zq plane, due to bending, 
the angle of inclination being 0 cos $ + ip sin $. 

The aerodynamic lift t, acting normal to the relative fluid 
velocity, and containing a term to account for the apparent camber 
due to rotation of the section about the x-axis, the drag d acting 
parallel to the relative velocity, and the moment ^ c y 2 acting about 

mid-chord, are given by 

£ = 2 p ^ 2 CC£(a) + |pVC 2 [ 4 - fi(e cos $ + ip sin $) ] 

a = ipv2cc d («> (6) 

M =/2 ' i ( ' v2c2c m < «> 


where 

V 2 = V 2 + V 2 . 
n t 
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I 


The blade chord is C and C^, C d and C m are experimentally deter- 
mined two-dimensional force coefficients. The apparent-camber 
term arises due to the linear variation of quasi-steady downwash 
along the blade chord when the blade section is rotated about mid- 
chord. The aerodynamic forces appearing in Eqs. (4) are related 
to those in Eqs. (6) by the following expressions: 


F = - t cos a - d sin a 
v 


F = £ sin a - d cos a 
w 


(7) 


M, = M . - Z l 

<j> c /2 a • 


The next step is to perform a consistent linearization of 
Eqs. (7) by expanding the various functions appearing in those 
equations in Taylor series and discarding second-order quantities 
in v, w, <t> , 0 or i(i. The expansion is made about the nominal inci- 
dence ao and nominal relative speed V 0 , where 


a 0 = $ + C , 


Vj- sin a - w. 

f s 1 

fix + V_ cos a sin fit ' 

f s J 

\ / 2. 

Vo = [ (v f sin a - w. ) 2 + (fix + V,. cos a sin fit) 2 ] ' ; 

r si ns 

as obtained when all dependent variables vanish in the expressions 
for V and a. Specifically, the expansions of a and V yield 

a = ciq + Aa + . . . 


£ = tan 


= a o + <f) + [ ftx + V, 

vS ‘ 


- w sin $ 


cos a sin fit] { [v + Z <f> 
s a 

- fi ( Z - JL ) (0 cos $ + sin $)]cos 
a z 

+ cos c* s cos fit (6 cos $ + ijj sin 4>) 


$ 

} 
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1 • • • 

- — [V^ sin a - w.]{w cos $ + [v + Z 41 
y2 1 s 1 a 

0 

- fi(Z - £ ) (0 cos $ + \p sin $) ] sin $ } + ... 

cl Z 


V = V 0 + AV + . . . 

= v ° + *- V f sin “s “ w i^ { tv + Z a 4 

- fi(Z - £ ) (e cos $ + if) sin $)]cos $ 

cL Z 

- w sin $ + cos a g cos fit(0 cos $ + ip sin $) } 


+ [fix + Vj. cos a sin fit] {w cos $ + [v + Z d> 
Vq r s 3 


- fi (Z - £ ) (0 cos $ + \p sin $)]sin $} + ... 


Each term in Eqs. ( 7 ) is then expanded, using these expressions 
for a and V. For example, the first term in is expanded by 

writing £ cos a as follows: 


£ cos a = {^pV 2 CC» (a) + 7pVC 2 [<J> - fi ( 0 cos 4 > + \p sin $) ] }cos a 


= {^ptV 2 + 2 V o AV + ...]C[C,(a 0 ) + 


dC, 


'£ ya ° J T dT 


A cx * 4 “ • • • ] 

ot = a 0 


+ ? P tV 0 + AV + . . . ] C 2 [ <£ - fi (0 cos $ + ijj sin <S>) ] } [cos a ( 


— sin ag Aa + . . . ] 


Thus , on discarding higher-order terms and the steady state term 
ipVgCC, (ag)cosao , which is not of concern in a stability analysis, 

£ cos a is given by 
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i dc e r 

Z cos a = 2 pC ^~ da a = ao cos “° ~ c £( a o) sin 

+ [fix + V, cos a sin fit] { [v + Z £ - fi(Z - Z ) (0 cos $ 

IS a a Z 

+ \p sin $)]cos $ 

- w sin $ + V* cos a cos fit(0 cos $ + i|> sin $} 

f s 

- [V.p sin a - w. ] {w cos $ + [v + Z £ - fi(Z = - £_) (0 cos $ 

I S X a. a. Z 

+ ip sin $) ] sin $ }| 

+ p CC » ( a o ) COS ao ■( [V- sin a - w. ] { [v + Z — fi (Z - Z ) (6 cos $ 

x- lx SI 3. 3 Z 

+ iJj sin $ ) ] cos $ 

- w sin $ + Vj cos a g cos fit(0 cos $ + 4 sin $) } 

+ [fix + V, cos a sin fit] {w cos $ + [v + Z £ - fi(Z - Z )(0 cos $ 

x S 3 3 Z 

+ 4/ sin $ ) ] sin $ } j 

+ ^pVoC^ cos ot o 1 ~ 0(0 cos $ + \p sin $) ] 

It can be seen at this point that, by systematically and con- 
sistently expanding all terms in the manner outlined above, the 
effects of compressibility, stall and reversed flow have been 

retained in the formulation. The coefficients C e , C, and C and 

Z a m 

their derivatives with respect to incidence, all evaluated at the 
nominal incidence ao(x,t), appear in the formulations as peri- 
odically varying factors in the expressions for the coefficients 
in the equations of motion. Thus, provided the correct and 
complete variations of those coefficients with nominal incidence 
and nominal Mach number are retained, the effects of compressi- 
bility, stall and reversed flow on the linear stability of the 
system are properly taken into account. 

Once Eqs . (7) have been expanded as outlined above, those re- 

lations are substituted in Eqs. (4) , which are in turn substituted 
in Eqs. (2) . With the blade motions expressed in terms of the 
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generalized coordinates, through Eqs. (1), the perturbation equa- 
tions of motion are obtained. Specifically, those equations are: 


N 


J k + 5 k e k'i l tGl “ ttUn + Bl “ <t,C n 1 = ° 

k — 1,2,... 


( 8 ) 


where 


R 


G to (t) - 


{ 4 K> + 4*’ (x)l ‘m (x ' t) 

+ a ! 9 (x) n (x,t) } dx 


R 


H kn (t) * 




{Ai k) (x) A,_ (x, t) + A^ k) (x)A (x,t) 


vn 


w wn 


+ Af k) (x)A_(x,t) + A| k) (x)x >I>n U,t) dx 


<j> ' ' 4>n 


> " ipn 


u.„(x,t) = 5 c «„ A, (n) 


vn 


D 4> 


+ | V 0 C [cos $A^ + sin $A^ Il/ ] {° D t v f cos ? cos a s cos 


(n) 


fit 


- fl(Z - X z ) cos a 0 ] + 2a L [fi(Z & - Z % ) sin a 0 


- Vjr sin £ cos a cos fit] + 5 C fi cos aQ} 

£ S 


V vn 


(x,t) = - 2m.fi e sin $A* n) - fpV 0 C 2 cos a 0 A 


(n) 


+ § V 0 C{[a D cos a 0 “ 2 ® L sin a o3 t z a A J n) + A v * J 


- [o n sin ao + 2 o L cos a<)3 A w ^ ^ } 
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in which. 


*D = 


dC. 


V«o> - d T 


a=ao 


sin ap " 


dC i 

c d(“°) + dct 


a=a 0 


cos ao 


= 2 jc^(a 0 ) cos a 0 + C d (a 0 ) sin ao] 


wn 


(x,t) = f VU t l A, 


(n) 


2 o ’h 9 
(n) 


(n) 


+ f Vo c A^ n; COS $ + a;“ # sin $ 


Vj; cos a - COS ? 
X s 


- f2(Z a - l z ) cos a 0 + T d 


ao 


V_ cos a sin % cos Sit 
f s 


- Si(Z - £ ) sin ao I - ^ C Si sm ao 

cl 2 


TT 

o * 2 


i n) + 1 pv 0 c 2 sin «„A< n) 


i m (x,t| = - 2 iafle cos $ +4 

Y l cos ao + Y d s; *- n ao 
- y l sin ao + Y d cos a 0 


+ ^ Vq C 


(2 &‘ n) + A‘ n) ) 
a 0 v 


, (n) 
w 


in which 
Y L = \ C l M 


dc d I 1 r , , dC t 

- ~s — ~ 1 cos a 0 + |C d (a 0 ) + da 

da I a=ao I J a 


a=ao 


sm ao 


D = 2 jc £ (a 0 ) sin a 0 “ C d (a 0 ) cos a^j 


cos Sit 
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* § V o c 6lR J n> 


+ f V 0 c 




Vj. cos a cos fit | 6 1 cos 5 + 26 sin 

X S 


- fi(Z & - £ z ) JjS 1 cos ao + 26 sin aoj 


+ 5 fi CZ V I cos 4> + A sin $ 
2 a|| 0 ip 


■] 


(n) n 


X *n (x ' t} = 2I 0 fl sin * % " 4 P V ° C Z a A <|> 


(n) 


+ § V 0 C 


6 1 cos ao + 26 sin ao 


0 Z A 
u |_ 


(n) (n) 

4> v 


- *i 


6 1 sin ao + 26 cos ao 





in which 


6 = C ( a o ) - % a C^(a 0 ) 


dC 


6 1 = C 


m 


da 


dC / 

_ z i 

a=ag a da 


a=a 0 


X (x, t) = 
ipn 


2l o fiA^ n ^ + 2mfie 


A^ n ^sin $ + A^ n ^ cos $ 
v w 
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DEVELOPMENT OF THE BASIC METHOD 


The specific problem at hand here is to establish the means 
for determining the stability of a set of linear equations with 
periodic coefficients , such as those representing a helicopter 
rotor in forward flight, Eqs. (8). Those equations are put in a 
somewhat more convenient form by defining an independent variable z, 

z = |ftt , 


and changing the notation of the coefficients, according to 


a mn ft H mn 


mm ft 


mn 


to 2 - G 
m mm 


ft 2 mn , m^n 


> m,n = 1,2,...,N. 


The equations to be analyzed then become 


a 2 ? 


m 


N 


. 2 + 1 

dz z n=i 


d^n 

a mn dz + ^mn *°n 


= 0 


(9) 


Ttl — 


where a and b are periodic functions : 
mn mn ^ 


+ ") = a mri (z) 
mn mn 


b (z + it) = b (z) 
mn mn 


As was noted in the introduction, equations of the form of 
Eqs. (9) have received considerable attention over the past cen- 
tury. Although attempts at general solutions have met with little 
success, the form of the solution in the general case has been 
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developed. The theory of Floquet (see Reference 6) yields that 


C n (z> = e la)Z <j> n ( z) 


where <t , n (z) is periodic, with a period it, and iw is a complex con- 
stant. Since the differential system (Eqs. (9)) is of order 2N, 
there are 2N values of ioo defining the solution for a given case. 

If any one of these characteristic values has a positive real part, 
the system is unstable. 

In order to secure the theoretical basis for the solution and 
to avoid numerical difficulties , it is necessary to first operate on 
Eqs. (9) to obtain two related sets of equations. It is convenient 
for this purpose to use matrix notation. A column m-'-'-rix X and 
square matrices A and B can be defined as 



r ?l' 


'ai l 

a i 2 • • • 

a iN 


r bn • 

b 12 

**• b lN 


£ 2 


a 2 i 




b 2 i 



X = 

• 

• 

A = 

• 

• 



B = 

• 

• 




• 


* 

3 N2 



• 

[ b N 1 

b 

N 2 

• • * b NN ; 


so that Eqs. (9) can be written 


d 2 X 

dz 2 


+ 



+ BX 


0 


( 10 ) 


where the derivative of a matrix is the matrix formed by replacing 
each element by its derivative. 

The first step in obtaining the two related sets of equations 
is to differentiate Eq. (10) , yielding 



(ID 
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If Eq. (10) is solved for d 2 X/dz 2 and the result substituted 
in Eq. (11) , it is found that 


d 3 X 

dz 3 


+ C 


dX 

dz 


+ DX 


0 


( 12 ) 


where 


C 



B-A 2 


D 


dB 

dz 


AB . 


If, now, Eq. (12) is differentiated and the second derivative 
of X eliminated as before, it is found that 


where 


d 4 X 
dz 4 


+ E 


dX 

dz 


+ FX 


0 


E = 


dC 

dz 


+ D - CA 



CB . 


(13) 


It can be shown that a set of functions is a solution of the 
original equation (Eq. (10)), if and only if it is a solution to 
both of the derived equations, Eq. (12) and Eq. (13) . As a result, 
the solutions of Eq. (10) can be found by solving Eqs. (12) and (13) 
and identifying those solutions common to the latter two equations. 
The proof is straightforward, and is given in Appendix A. 

Consider, first, the solution to Eq. (12) . It is convenient 

at this point to abandon the matrix notation. Thus, if c and 
c mn 

d mn denote the elements of matrices C and D, respectively, Eq. (12) 
gives that 


d 3 c N 

51 + l 

dz 3 n= l 


r 


mn 


d £ 

dz + mn 


n 


0 


(14) 


m = 1,2, ...,N . 
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From the theory of Floquet, the solution of Eqs. (14) must be 
of the form 


V z ’ ■ el "V z > 



p mk e i(2k+ “ )z 


m *“ 1 , 2 , * • • , N f 


upon expansion of <|> m in a complex Fourier series. Also, the peri- 
odic coefficients in Eqs. (14) can be expanded in Fourier series: 


c = T £ , e 

mn i L mnk 
k=-°° 


2ikz 


l n e 


2ikz 


mn , u 'mnk 
k=-<» 


}■ m,n = 1,2 , . . . ,N, 


If the Fourier representations for the solution and for the 
coefficients are substituted in Eqs. (14) and the coefficients of 
2. ik z 

e are grouped and set equal to zero, a set of linear recursion 

relations in the unknown coefficients p m ^ is obtained. Specifi- 
cally, it is found that 


0 = [i ( 2n + a>) ] 3 p rn 


+ 


N 

I 


S=1 



n rsk 


+ 


i (2n - 2k + w) 




(15) 


n — ..., 2, 1,0,1,2,...; 
r = 1 ,2 , . . . ,N 


20 



These relations constitute an infinite set of linear algebraic 
equations in the unknown Fourier coefficients p^. For this set 

of equations to have a nontrivial solution, the associated infinite 
determinant A (to) must vanish (a discussion of infinite determinants 
is given in Reference 6) . The requirement that A vanish is the con- 
dition which determines to, and hence the stability of the system. 

In order that A be meaningfully defined, it is necessary to 
divide each of Eqs. (15) through by the coefficient of P rn * With 

the unknowns then appropriately ordered, the diagonal elements of 
A are all unity and the off-diagonal elements all vanish in the 
limit as a row or column index of the determinant tends to either 
positive or negative infinity. Specifically, the elements of 

A ( to) ( y , v = 0, ±1, ±2...) are given by 


cr„ . , = a (n,k;co) 

Nn+r, Nk+s rs 

n,k . • . , — 2, — 1 , 0 , 1 , 2, ...7 

r ,s — 1,2,...,N, 


whe re 


a rs (n,k;w) = 


[i(2k + »)] 3 S nk 5 rs + i(2k + ») S rs-n _ k * " rs , n - k 

[i(2n + w) ] 3 + i(2n + to) E + n 

rr 0 rro 


in which 


<5 i j = 0 , i / j ; 


= 1 , i = j. 


Note that the diagonal elements a are all unity. The con- 
struction of the determinant can perhaps be best visualized as 
made up of a collection of subarrays which are NxN in size. The 
location of any element within a subarray is determined from 
indices r and s. The location of each subarray is specified 
through indices n and k. The general arrangement of the a's in 
A is illustrated in Figure 2 for the case N=2. 

The value of A is obtained, for a given to, by evaluating the 
finite determinant formed by ranging n and k from, say, -L to L. 
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1 

“12 (-1,-1) 

“1 i (“1,0) 

a 2 i (-1,-1) 

1 

“2 1 (-1,0) 

“11 (.0,-1) 

“12 (0,-1) 

1 

“2 1 (0,-1) 

“22 (0,-1) 

“21 (0,0) 

“1 1 (1,-1) 

“12 (1>-1) 

“11 (1,0) 

«2l(l/-l) 

“22 (1,-1) 

“2 1 d,0) 


“ 12 (- 1 , 0 ) 

“ 11 (- 1 , 1 ) 

“ 12 (“ 1 , 1 ) 

“ 22 (“ 1 , 0 ) 

“ 2 1 (- 1 , 1 ) 

“22 (- 1 ,D 

“12 ( 0 , 0 ) 

“11 ( 0 , 1 ) 

“ 12 ( 0 , 1 ) 

1 

“21 ( 0 , 1 ) 

“22 ( 0 , 1 ) 

“ 12 ( 1 , 0 ) 

1 

“ 12 ( 1 , 1 ) 

“ 22 ( 1 , 0 ) 

“2 1 ( 1 , 1 ) 

1 


Figure 2 . The arrangement of the determinant 
elements for the case N = 2 


Successively larger values for L are then taken until the limit be- 
comes apparent (see Reference 6) . 

The expression 

A (to) = 0 (16) 


by itself is not a particularly useful relation for determining 
w, because of the limiting process involved in evaluating infinite 
determinants. It will be shown, however, that A(w) in fact con- 
stitutes a combined series-product expansion in w of a finite sum 
of trigonometric functions. With a (oj) expressed in the latter 
form, Eq. (16) can be solved explicitly for w. 

It should be noted here that Hill obtained such a relation 
for the second-order differential equation which bears his name 
(see Reference 6) . The development which follows is a generaliza- 
tion of that result. 

Two properties of the function A(io) must first be established. 
Specifically, it is asserted that 

1. A (to) is an analytic function of w, except for simple, 
identifiable poles; 

2. A (<d) is periodic in to with a period of 2. 

That A (to) is analytic can be concluded by noting that A is an 
absolutely convergent determinant — i.e. , the product of the diag- 
onal elements converges absolutely (in this case to unity) and the 
double sum of the off-diagonal elements converges absolutely. 
Uniform convergence and analyticity can then be established (see 
Reference 6) . 

Note that if the original equations (Eqs. (9)), rather than 
derived equations (Eqs. (14)), had been used to generate A, that 
determinant would not have been absolutely convergent since the 
double sum of off-diagonal elements generated from the original 
equations is only conditionally convergent. What follows hinges 
on the analyticity of A , hence the necessity for working with the 
differentiated sets of equations. 

Clearly, too, the only singularities in A are simple poles at 
the points 


b) 


2n - iX . 


r 3 


n = 0 , ±1, ±2 , . 


r = 1,2, . . .N; 
j = 1/2,3; 


23 



where the X 1 s are the 3N roots of the N cubic equations 


+ X? + n — 0 , r — 1,2 , . . . ,N. 
rrg rrg 


The periodicity of A is shown by substituting to + 2 for to in 
the expression for a , whereupon it is found that 

jT S 


a (n ,k ; co + 2) = a _ (n + l,k + l;to) 

J- S i s • 


Thus, in the limit, the value of A is unchanged if to is replaced 
by a) + 2 : 


A ( to + 2) = A (to) . 


Now, consider the function D(to), defined by 


N 


D(to) = A ( to ) + y y f . cot[S(to + ix .)] 

<-> . T I / T I 


r= l ] = l 


r: 


rj 


(17) 


where the f r j ' s are constants. Observe that D(iu) is an analytic 
function and that 

D(oj + 2) = D(to). 


Further, note that the term 


f . cot [? (to + ix . ) ] 


has simple poles at to = - 2n - iX r j for n = 0, ± 1, ± 2, 


and has 


no other singularities, 
is properly chosen, D(co) will have no poles. 


Thus, if the value for each constant f . 

r D 

Let the f .'s be so 


rj 


chosen, making D(to) analytic throughout the finite part of the <o- 
plane. But D is clearly bounded at infinity; A (to) has a limit of 
one and the cotangents have limits of ±i for Im ( to) ->- ±o ° . Therefore, 
by Liouville's theorem, D(u>) is simply some constant, say c: 
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D ( a)) = c = A ( 03 ) 


(18) 


N 


+ l l f r . cot[|(to + iA rj >] . 


r= l j = i 


r: 


It only remains to determine the values of the f r ^ 1 s and of 

c. This is facilitated by first considering the limits of D(to) 
with infinite to: 


N 3 

lim DU) = c = 1 - i l l f. f 

Im(to)-*- + °° r=i j = i J 


lim 
Im( to) ■ 


N 3 

DU) = c = 1 + i £ l f. 
0 r= i j = i J 


Clearly , then, c=l and 


N d 

I I 

r=i j=i 


f rJ - 0 


Replacing c by unity and solving for A in Eq. (18) gives that 


N 


A ( to) = 1 - l l f r j cot[|(to + iX r _. ) ] . 


r= l 3 = 1 


r 3 


Now, let 3N - 1 arbitrary (but finite) values of to, say toi, 
t 0 £ , . . . , W 3 N _ 2 _ , assigned in the above equation. The resulting 

3N - 1 equations, together with the one obtained for infinite to, 
provide sufficient relations to solve for the f^'s. More speci- 
fically, those constants are the solution of 

N 3 

£ £ f ■ cot [tU-^. + iX . ) ] = 1 — A( ton,) t ^ = 1,2,..., 3N — 1 } 

r=i j = i r3 Z K 303 K 


N 3 

l l f 

r=i j=i 


r 3 


= 0 . 


(19) 
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With the f r j 1 s known, the determinantal equation, A (to) = 0, 
can be replaced by 


N 3 

1 - I l f , cot[^ (u + iX .) ] = 0 . (20) 

r=l j=i r3 rD 


A polynomial of degree 3N in e 17rw can be readily derived from 
Eq. (20) . The 3N roots of that polynomial determine the charac- 
teristic values for Eqs. (12). 

In the same manner as is outlined above for Eqs. (12) , the 
4N characteristic values for Eqs. (13) can be obtained. The 2N 
values common to the two sets are those of the original equations, 
Eqs . ( 9) . 


OUTLINE OF BASIC COMPUTER PROGRAM 


In order to implement the method derived for analyzing sta- 
bility of dynamic systems with periodic parameters, a basic digital 
computer program was prepared for the analysis of systems with 
three degrees of freedom. The program accepts the pertinent data 
for a given system in the form of Fourier coefficients of the 
periodic coefficients in the equations of motion (Eqs. (9) , with 
N=3) , and then proceeds to calculate the six characteristic 
values of the system by the method derived in the previous section. 
For specific applications, a small subroutine can be prepared, if 
necessary, to generate the Fourier coefficients used as input by 
the basic computer program. 

The overall flow of information guided by the formulations 
of the basic program are outlined below. The numbers in the out- 
line correspond to those in the blocks on the flow diagram of 
Figure 3. The outline includes reference to the pertinent equa- 
tions of the previous section. The relations guiding the more 
routine procedures, such as extraction of the roots of polynomials 
and solving of linear algebraic equations, have been omitted, but 
are embodied in the computer program, a listing of which is given 
in Appendix B. 

1. Fourier Expansion of Coefficients 

Given the coefficients of the equations of motion for a linear 
dynamic system with three degrees of freedom and periodically 
varying parameters with normal modes used as generalized co- 
ordinates, such as Eqs. (9), the coefficients are expanded 
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Figure 3. Procedure for calculating characteristic values of a 
dynamic system with three degrees of freedom and 
periodically varying parameters 
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in Fourier series. The specific output is a set of Fourier 
coefficients for the periodic coefficients and fc> mn appear- 

ing in'Eqs. (9) with N=3; m,n=l,2,3. This subroutine is 
associated with the particular dynamic system being analyzed, 
and can be described as preparing the input to the main 
stability-analysis program. 

2. Matrix Multiplication and Fourier Expansion of 
Coefficients 

Given the Fourier coefficients of the elements of the 3x3 
matrices A and B appearing in Eq. (10) , the Fourier coeffici- 
ents of the elements of matrices C and D, appearing in Eq. 

(12), and of E and F, appearing in Eq. (13), are calculated. 

In this way, the coefficients of the higher-order, or differ- 
entiated, systems of equations are related to the original 
system of equations, Eqs. (9). 

3. Polynomial Factorization 

The three cubic polynomials associated with the determinant 
for the ninth-order system are factored to obtain the nine 
singularities of that determinant (subroutine 3a) . The three 
quartic polynomials associated with the twelfth-order system 
are factored to obtain the twelve singularities of that deter- 
minant (subroutine 3b) . The points w k at which the infinite 

determinants are to be evaluated are also selected, by speci- 
fying each of them to be a set distance <S from one of the 
singular points of the determinant. The accuracy of the solu- 
tion was found to be quite sensitive to the value of 6, it 
generally being necessary to make 6 as small as possible, 
without sacrificing numerical accuracy elsewhere in the pro- 
gram, in order to get satisfactory results. The reason for 
this behavior is not known, but is presumably connected in 
some way to the accuracy with which the infinite determinants 
are calculated, this being the most difficult, and hence 
least accurate, of the tasks performed by the program. A 

— 4 

value for 6 of 10 was used for most of the calculations 
reported here. 

4. Determinant Element and Determinant Evaluation 

For each value of oj^, the determinant elements are computed 

and the infinite determinant is evaluated. It was anticipated 
that this step would be the most time-consuming one in the 
program, and so considerable care was taken to employ the 
most economical means for determinant evaluation. A process 
of tri angulari z at ion was selected for evaluating finite 
determinants. 


28 


II 



as the 


The limiting value of a determinant, for a given 

number of elements increases without bound, is estimated as 
follows. Three determinant values, for three successively 
larger numbers of rows and columns, are first obtained. The 
determinant is then assumed to vary inversely as some unknown 
power of the number of rows and columns. The assumed form 
has three unknown constants , so the three determinant values , 
corresponding to three different determinant sizes, provide 
sufficient information to calculate the constants and hence 
an estimate of the limit as the number of elements becomes 
infinite. It has been found that three determinants, re- 
spectively 33x33, 39x39 and 45x45 in size, generally are 
sufficient to provide three-place accuracy in the result. 

5. Set-up and Solution of the Linear Algebraic Equations 


For the ninth-order system, the eight values of and the 
eight values of A (to^.) are used to compute the coefficients 

and inhomogeneous terms of Eqs. (19) . The nine equations 

are then solved for the nine constants f .. A similar pro- 

r D 

cedure is followed to obtain the twelve constants for the 
twelfth-order system. The equations are solved by triangu- 
larization. 


6. Determination of the Characteristic Polynomials of the 
Higher-Order Systems 

For the ninth-order system, the coefficients of the poly- 
nomial of the ninth degree in e ilTW are derived from the deter- 
minantal equation, Eq. (20) (subroutine 6a) . Similarly the 
coefficients of the characteristic polynomial of the twelfth- 
order system are obtained from the appropriate determinantal 
equation (subroutine 6b) . 

7. Evaluation of the Common Polynomial Factor 

A.t this point, one could presumably compute the nine roots of 
the ninth-degree polynomial, the twelve roots of the twelfth- 
degree polynomial and identify the roots common to the two 
polynomials as being the characteristic values of the original 
system. However, this procedure requires evaluation of fifteen 
extraneous roots, followed by possible difficulties in identi- 
fying which roots are indeed common ones . 

It was found, though, that the coefficients of the character^ 
is tic equation for the original system, which is actually a 
polynomial factor common to the two higher-degree polynomials. 
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can be obtained in terms of the coefficients of these higher- 
degree polynomials. The steps taken to derive the necessary 
expressions are outlined in Appendix C. Subroutine 7 imple- 
ments those expressions, yielding the coefficients of the 
sixth-degree polynomial characterizing the original system. 

8. Polynomial Factorization 

A standard library subroutine is used to obtain the six roots 
of the characteristic polynomial. The logarithm of each root 
is then evaluated (recall that the polynomial is formed of 

powers of e lira) ) to obtain the six characteristic values, and 
hence the stability, of the system. 

During the check-out of the computer program, it was found 
necessary to extract the roots of the higher-degree poly- 
nomials. Since little additional running time is consumed 
by those calculations, determination of the roots and charac- 
teristic values for the ninth and twelfth-order systems has 
been retained in the program for purposes of comparison. 
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APPLICATION OF THE METHOD 


The computer program implementing the general method was first 
checked out by means of a test case, the characteristic values of 
which could be derived in advance. The program was used next to 
analyze a rotor system with two degrees of freedom for which solu- 
tions by direct time integration had been obtained previously. 
Lastly, calculations were conducted for a model rotor system having 
three degrees of freedom, for which flutter boundaries had been 
obtained experimentally. These applications are discussed in detail 
below. 


Test Case Calculations 


It was found necessary in the course of the check-out of the 
computer program to have available a test case for which the char- 
acteristic values of the original sixth-order system as well as of 
the derived ninth-order system were known in advance. Such a case 
was generated by appropriately transforming Mathieu's equation. 

The sixth-order system was then made up from three of these equa- 
tions, treated as a coupled system, the program not being able to 
distinguish whether the equations are coupled or not. 

The derivation of the test case was carried out as follows. 
Consider Mathieu's equation, which is generally written in the 
form 


+ (a - 2q cos 2z)u = 0 (21) 

dz 2 


where a and q are specified constants. The characteristic values 
±y of Mathieu's equation are a function of a and q (see Reference 


5) . 


Now, let a new dependent variable y be defined by 

- |f(z) 

y = ue 


where f(z) is periodic with a period n , but otherwise may be re- 
garded as arbitrary. If u is substituted in Eq. (21) in terms of 
y and f the following differential equation for y is obtained: 
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§hL + A §2. + By = 0 
dz 2 dz 


( 22 ) 


where 


A = | f ( z ) 


B = a - 2q cos 2z + f 2 + ~ — 

36 6 dz 


Suppose, now, f is given by 


oo 

f (z) = y (A cos 2nz + B sin 2nz) . 
u „ n n 

n=o 


It then follows that the characteristic values of y, as a solution 
of Eq. (22) , must be p - A 0 /6 and - y - A 0 /6. Thus, we have derived 

a second-order differential equation with periodic coefficients, of 
a quite general form, for which the characteristic values are known. 

Further, suppose a third-order equation is derived by the 
method described previously which eliminates the second derivative, 
namely 


+ C ^ + Dy = 0 (23) 

dz 3 dz 


where 


C 

D 


B + 


dA 

dz 


- A 2 


dB 

dz 


AB 


It is possible, through further transformation and utilization of 
the relations given on p. 134 of Reference 5, to show that the 
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additional characteristic value of Eq. (23) is the negative of the 
sum of the other two values, i.e., Ao/3. 

The three equations of the test case were all derived from the 
Mathieu's equation having a = 1, q = 0.32. For those values of a 
and q, it follows that y = .15813 + i (see p. 105 of Reference 5) . 
Only the first three terms of the series for f(z) were retained 
for each equation of the test case. The values of A 0 , Aj and A 2 

chosen and the associated characteristic values for the three de- 
grees of freedom are as follows: 


A 0 

Ai 

a 2 

A 0 

y ' ”6" 

A 0 

- y - — 

Ao/3 

- 4.5 

2.0 

0.6 

0.90813 + i 

0.59187 - i 

- 1.5 

- 0.6 

2.0 

- 0.6 

0.25813 + i 

- 0.05813 - i 

- 0.2 

3.0 

- 2.0 

0.6 

- 0.34187 + i 

- 0.65813 - i 

1.0 


The Fourier coefficients of the functions A and B appearing 
in Eq. (22) were calculated for the three degrees of freedom and 
were supplied as inputs to the computer program. The character- 
istic values which were calculated for the ninth and twelfth order 
systems are tabulated below, together with the independently 
derived exact (anticipated) values. The quantities A and A T 

R, Jl 

given are related to u> by 


iw = X R + iAj . 
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CHARACTERISTIC VALUES OF THE 9TH ORDER SYSTEM 


Calculated 

Anticipated 

X R 

X I 

X R 

X I 

0.90813 

1.0 

0.90813 

1.0 

0.59185 

- 1.0 

0.59187 

- 1.0 

0.25814 

1.0 

0.25813 

1.0 

- 0.05816 

O 

* 

i — I 
1 

- 0.05813 

- 1.0 

- 0.34184 

1.0 

- 0.34187 

1.0 

- 0.65813 

- 1.0 

- 0.65813 

- 1.0 

- 1.50000 

o 

• 

o 

- 1.50000 


- 0.19998 

o 

• 

o 

- 0.20000 

0.0 

1.00000 

0.0 

1.00000 

0.0 


CHARACTERISTIC VALUES OF THE 12TH ORDER SYSTEM 


Calculated 

Anticipated 

X R 

X I 

X R 

X I 

0.90803 

1.0 

0.90813 

1.0 

0.59187 

- 1.0 

0.59187 

- 1.0 

0.25802 

1.0 

0.25813 

1.0 

- 0.05960 

- 1.0 

- 0.05813 

- l.o 

- 0.34188 

1.0 

- 0.34187 

1.0 

- 0.65813 

- 1.0 

- 0.65813 

- 1.0 

0.97699 

1.0 

- 

- 

0.02369 

- 1.0 

- 

- 

- 0.10970 

1.0 

- 

- 

- 0.41425 

0.0 


- 

0.21431 

0.0 

- 

- 

- 1.38935 

- 1.0 

- 

- 
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As can be seen from the tabulated results, the characteristic 
values common to the ninth and twelfth order systems are those of 
the sixth-order system. Note, too, that the program for the most 
part predicts the characteristic values for this case to four 
significant figures . 

The characteristic values were also calculated from the sixth- 
degree polynomial which was derived from the two higher-degree 
polynomials. The resulting roots were unacceptably different from 
the correct values, apparently because the polynomial as so derived 
is sensitive to slight inaccuracies in the coefficients of the 
higher-degree polynomials for this case. Thus, it may be necessary 
to rely on direct comparison of ninth-order and twelfth-order solu- 
tions to determine the correct sixth-order solutions in some in- 
stances . 


Comparison With Direct Time Integration 


The analysis of stability by direct time integration on a digi- 
tal computer of a rigid rotor blade with flapping and lead- lag 
hinges is reported in Reference 12. The nonlinear representations 
of both inertial and aerodynamic forces are utilized in the equa- 
tions of motion. The basic blade for whxch numerical results were 
obtained had zero offset of the flapping hinge and 0.05R offset of 
the lead-lag hinge. The blade had a constant chord except for a 
cut-out from the axis of rotation to 0.2 R. In the calculations 
reported in Reference 12, the rotor was unloaded and at zero shaft 
angle. Further details can be found in Reference 12. 

For two of the cases analyzed in Reference 12, the variations 
of flapping and lead-lag angles with time are presented. These 
cases had values of advance ratio y of 0.6 and 1.4, respectively. 

The mass constant y 1 for both cases was 1.6 (y 1 = 0 CR 4 /I h , where 

1^ is the mass moment of inertia about the lead-lag hinge) . The 

case for the lower advance ratio is reported to be very stable and 
the case with y = 1.4 is indicated to be stable but near a boundary 
of neutral stability. The time histories of the blade displacements, 
from Reference 12, are reproduced in Figures 4 and 5. 

The appropriate parameter values for these two cases were in- 
serted in the linearized equations of motion of a rotor blade 
developed previously in this report. Series representations of the 
coefficients in the equations, including the first eleven harmonics 
of rotor rotational speed, were generated and supplied to the main 
stability-analysis program. The characteristic values calculated 
are as follows: 
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DEGREE OF 
FREEDOM 

CHARACTERISTIC VALUES 

U = 0.6 

y = 1.4 

X R 

X I 

a r 

X I 

Flapping 

- 1.13400 

- 1.13400 

0.60206 
- 0.60206 

- 0.73379 

- 2.56010 

- 1.0 
1.0 

Lead-Lag 

- 0.00353 

- 0.00353 

0.56195 
- 0.56195 

- 0.00534 

- 0.00534 

0.56206 
- 0.56206 


Examination of these results indicates qualitative agreement 
with the results of the analysis by direct time integration but 
there is evidence of some differences in quantity. At y = 0.6, the 

_ 2 

flapping motion should damp by a factor e = 0.135 in 4/(-A R ) 

= 3.527 radians of azimuth change, by the result obtained here, 
whereas Figure 4 indicates a much more rapid decrease. On the 

— 2 

other hand, at y = 1.4, the factor e should apply to the flapping 
motion for a change of 4 /(-a t J = 5.452 radians, but Figure 5 indi- 

.K 

cates that the flapping motion is considerably less stable than 
that. 


From the qualitative viewpoint, the comparison is more favor- 
able. The decrease in stability of the flapping degree of freedom 
with increasing y is in evidence in the result obtained here, 
since X R is less negative at y = 1.4 than at y = 0.6. The sta- 
bility of a given system is determined, of course, by the least 
negative, or most positive, value of X_. Also, in agreement with 

the indications of Figures 4 and 5, the lead-lag motion is only 

slightly damped, with a factor e 2 decrease occurring in 1,130 
radians at y = 0.6 and in 750 radians at y = 1.4. 

The quantitative differences in the predictions of the flap- 
ping motion can be attributed to the nonlinear effects which are 
included in the direct-integration solution, but which are absent 
from the formulations analyzed here. This can be seen as follows. 
With the rotor unloaded, as is the case here, the equations for 
rigid-body flapping and lead-lag motion become decoupled when 
linearized. Thus, any coupling of the motion which is detected 
can be attributed, in this case, to nonlinear dynamic-coupling 
effects. Since the motions plotted in Figures 4 and 5 were 
initiated by a disturbance in flapping, the considerable lead-lag 
motion must then all be due to nonlinear effects. Furthermore, 
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the flapping motion for y = 1.4 at the higher azimuth angles 
appears to have been excited by the persisting lead-lag motion, 
resulting in an apparent damping of the motion which is much less 
than would otherwise be the case. 


Comparison With Experimental Data 

Reference 10 reports the results of an experimental investi- 
gation of the flutter of a model helicopter rotor in forward 
flight. The model rotor had a single blade with a radius of four 
feet and with a flapping hinge through the axis of rotation. The 
blade had a constant chord of 3.5 inches and a root cutout of 6 
inches. Inertial and elastic properties of the blade are given in 
Reference 10. The blade was relatively stiff in torsion, but the 
control system was made flexible, so the primary contributions to 
blade motions derived from rigid-body pitch, flapping motions and 
deflections in the first flapwise bending mode. The main parame- 
ters of the study were advance ratio, control system stiffness 
and chordwise mass center location. 

The case selected for comparison had a ratio of first flap- 
wise bending frequency to^ to control frequency to 0Q (nonrotating) 

of 1.31 and a chordwise mass center location of 42.5% of chord aft 
of the leading edge, giving a value of 0.139 for equivalent mass 
center location as defined in Reference 10. This case was chosen 
because the data as plotted in Figure 8e of Reference 10 indicated 
there should be a relatively large change in stability with ad- 
vance ratio. Natural frequencies and mode shapes for the first 
three coupled modes of the blade were calculated for a value of 

the ratio to Q /ft = 4.37, and Fourier coefficients of the coeffi- 

e o 

cients of the equations of motion were calculated for values of 
advance ratio y of 0.0, 0.09,. 0.175, 0.24 and 0.30. The Fourier 
coefficients for each value of y were supplied to the main com- 
puter program and the system characteristic values computed. 

Subsequent to these calculations , it was determined through 
communications with one of the authors of Reference 10 that the 
data of Figure_8e of that reference were mislabelled. The symbols 

for values of to, /to„ of 1.31 and 0.63 were reversed. The experi- 
$ 1 0o 

mentally determined flutter boundary actually corresponding to the 

calculations performed, consisting of a plot of to. /ft versus y, 

0o 

taken from Figure 8e of Reference 10, is reproduced in Figure 6. 
The points at which characteristic values were calculated are 

indicated by asterisks on the line drawn at to. /ft = 4.37. 

0o 
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As can be seen from Figure 6 , the mislabelling of the data 
has led to a rather unsatisfactory comparison of theory with 
experiment. The experimental data shows very little change with y. 
Hence, for a comparison with this case, the characteristic values 

should have been determined at several different values of w /fl, 

e o 

for fixed y, rather than at fixed to. /fl for various values of y. 

o o 

Unfortunately, time limitations prevented carrying out more ex- 
tensive calculations. However, some information can still be 
derived from the calculations which were performed. 

The characteristic values obtained for each advance ratio 
are as follows: 


P 

= 0.0 

P 

= 0.09 

P 

= 0.175 

X R 

X I 

X R 

A I 

X R 

A I 

- 0.06273 

0.02028 

- 0.06271 

0.02800 

- 0.06266 

0.02871 

- 0.06273 

- 0.02028 

- 0.06271 

- 0.02800 

— 0 . 06266 

- 0.02871 

- 0.33052 

0.20163 

- 0.33055 

0.20185 

- 0.33068 

0.20267 

- 0.33052 

- 0.20163 

- 0.33055 

- 0.20185 

- 0.33068 

- 0.20267 

- 0.47015 

0.08874 

- 0.47016 

0.08954 

- 0.47012 

0.09152 

- 0.47015 

- 0.08874 

- 0.47016 

- 0.08954 

- 0.47012 

- 0.09152 


P = 

0.24 

P = 

0.30 

X R 

X I 

X R 

X I 

- 0.06268 

0.02957 

- 0.06283 

0.03060 

- 0.06268 

- 0.02957 

- 0.06283 

- 0.03060 

- 0.33120 

0.20366 

- 0.33227 

0.20489 

- 0.33120 

- 0.20366 

- 0.33227 

- 0.20489 

- 0.47010 

0.09342 

- 0.47013 

0.09510 

- 0.47010 

- 0.09342 

- 0.47013 

- 0.09510 
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natural frequency to Q to rotational speed J2 
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versus advance ratio u, for w. /to _ = 1.31 
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and x g /C = f 139 (from Reference 10). Points 

at which characteristic values were calculated 
are indicated by asterisks. 



The results of the calculations indicate, first, that the stability 
of the rotor for the control stiffness selected is essentially in- 
dependent of advance ratio over the range of y considered. The 
relative insensitivity to advance ratio changes is certainly in 
agreement with the plot of Figure 6. Secondly, it should be noted 
that the rotor is predicted to be very nearly unstable. The value 
for X of - 0.063 indicates that about ten rotor revolutions are 

R __ r\ 

required to damp the motion by a factor e . The limited calcula- 
tions which were performed, then, are at least in qualitative 
agreement with the experimental results. A definitive quantitative 
comparison must await more extensive calculations. 


CONCLUDING REMARKS 


A method has been developed for analyzing the aeroelastic 
stability of helicopter rotors in forward flight. The method 
employs formulations for calculating the characteristic values 
of the perturbation equations of motion, the latter being a 
coupled set of second-order, linear differential equations with 
periodic coefficients. The characteristic values, which are the 
zeros of an infinite determinant, are calculated from an equiva- 
lent analytic form for the infinite determinant consisting of a 
finite sum of trigonometric functions. 

A digital computer program was prepared which implements the 
method for three degrees of freedom. Calculations of charac- 
teristic values carried out for a test case demonstrated the 
practicality of the method, with anticipated results generally 
obtained to four significant figures. 

Calculations were carried out for comparison with results of 
a direct time integration on a digital computer of the equations 
for a rigid rotor blade with flapping and lead-lag hinges. The 
results were in qualitative agreement. Quantitative differences 
were attributable to the nonlinear effects which were included in 
the direct time integration. 

Lastly, limited calculations were performed for comparison 
with experimentally derived flutter boundaries for a model rotor 
blade with three degrees of freedom. The calculations indicate 
that the rotor is only marginally stable at the control stiffness 
selected and that the stability is relatively insensitive to ad- 
vance ratio, in agreement with the experimental results. 
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APPENDIX A 


Relationship Among the Solutions of 
Original and Differentiated Systems 


Let Xo be a solution of 

d 3 X 


o dx 0 

+ C + DX 0 = 0 

dz 3 dz 


as well as of 


From Eq. (A- 2) , 


d 4 X 0 dX 0 

+ E + FX 0 = 0 

dz **■ dz 


d 4 X 0 

d 2 A 2dB 

AB - — (A 2 ) - 

dz 

fdA _ 2 ] 

A 

dz 4 

jiz 2 dz 

, dz 

rA 


4 - 

P B - d ( AB) - 

dz 2 dz 

— + B - A 2 
^dz 

B 


dX 0 

dz 

Xn = 0 


But Eq. (A— 1) , differentiated once, gives 


d 4 X 0 


dz' 


— + B - A 2 
dz 


d 2 X 0 


dz 2 


+ dB _ d_ (A 2) + dB 
dz 2 dz dz dz 


dX 

dz 


d 2 B d 
=-=■ - — (AB) 

dz 2 dz 


X f 


( A— 1) 


(A- 2) 


(A- 3) 


0 


( A-4 ) 
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If (A-4) is substituted in (A-3) , a number of terms cancel and a 
common factor can be extracted, with the result that 


— + B - A 2 
dz 


d 2 X 0 dX 0 

+ A + BX 0 

dz 2 dz 


0 


Thus, X o must be a solution of the original equation, Eq. (10). 

It is, of course, also true that every solution of Eq. (10) is a 
solution of Eq. (A-l) and of Eq. (A-2) . Thus, Xq is a solution 

of Eq. (10) if and only if it is a solution of both Eq. (A-l) and 
Eq. (A-2) . 


44 



n n □ n n 


APPENDIX B 


Listing of Basic Computer Program — < = s£ez^ 

The program was coded in FORTRAN IV. A CDC 6400 digital 
computer was employed for all calculations. 


PROGRAM Q( INPUT, OUTPUT, T APE5* I NPUT , TARE6*0UTPUT ) 

DIMENSION L2(5> 

DIMENSION ER( 350 > ,EI ( 350 ) ,FR(350 > ,Fl <350 ) , FFR ( 150 ) , FF I < 150 ) 
DIMENSION DELT( 150) .ALPHA (ISO ) . BET A ( 150 ) , GAMMA ( 150 ) 

DIMENSION GR(5000), 31(5000) 

DIMENSION CI( 25 0),DR(250), 01(250), DuM(8l0 ) 

DIMENSION YR9(9), Y19(9), YR12(12), Yll2(12),BlGA(9) ,RIGAH(13). W(13) 
DIMENSION AK9(9).AK12(13),0YR9(11).DYI9(11),DYR12(12),DYU2(12) 

D I MENS I ON S I GMA ( 7 ) 

DIMENSION XI9(10).ET9(1Q),XI12<13).ET12(13),XI6(7),ETA6<7) 
DIMENSION R12TR(13) .R12TI (13) 

DIMENSION A UM (1, 1),AY(15),AR(125).AI(1 25 ),BR (125 ). 81(125), CR (250) 
DIMENSION PSI9(9),R06TR(6),R06TI(6),NDR9(9),NDI9(9),N DR 12(12) 
DIMENSION NDI12(12)»DETR<12),QETI(12),PIRNEW(12),PIR0LD<12> 
DIMENSION PI INEW( 12) , P I I OlO ( 12 ) , AUMR9 ( 8 , 14 ) , AUM I 9 ( 8 , 14 ) 

DIMENSION AUMR12(11,14) , AUM 1 12 < 11 , 14 ) 

DIMENSION RQLY(9),P5U2(13),=TA9(1Q)»ETA12<13) 

DIMENSION RR(2Q),RRI(2Q),ROOTM20),ROOTR< 20>,Z(20).Y<2Q).P3L¥(20) 
DIMENSION XR ( 12 ) »XI (12) 

DIMENSION CCH(150).CC(150),PSI6(7) 

DIMENSION RR12U3), RRU2(13) 

DIMENSION A9(9),B9(9),C9(9),AH(12),BH<12),CH<12) 

DIMENSION I CHG ( 15 ) 

THE GR AND GI ARRAYS ARE DETERMINANT ELEMENT STORAGE NOT NEEDED 
THROUGHOUT THE PROGRAM AND ARE AVAILABLE FOR FOR OTHER USES 3EF0RE 
AND AFTER THE DETERMINANT ROUTINE 

EQUI VALENCE ( ALPHA ( 1, 1 ) , AR ( 1 ) ) , ( BETA ( 1 , 1 ) , A I (1 ) ) , ( GAMMA ( 1,1 > ,3R(1) ) 
FQUIVALENUE ( UELT ( 1 , 1 ) , 31 ( 1 ) ) 

EQUI VALENCE (R12TR(1 ) , DU M < 1 )) . ( R 12 T I ( 1 ) , DUM ( 15 ) ) ( A Y ( 1 ) ,DUM(30 ) ) 

EQUI VALENCE (PS I 9( 1 ),DJM(50)).(RQ6TR(1),DUM(60)),(R06TI(1),DJM(70)) 
EQUI VALENCE(NDR9< 1) ,l)JM( 80 ) ) , (NDI9< 1 ) , DUM( 90 ) ) , (NDR12(1) , UUM( 100 ) > 
EQUIVALENCE (OYR9(l),DUM(H5)),(DYl9(l),DUM(130)),(RR(l),DU'K145)) 
EQUIVALENCE (DYRl2(l) , 0 U M ( 1 6 6 > ) , ( NDl 12( 1 > » DUM< 180 ) ) 

EQUI VALENCE <DETR( 1) ,IDJM(193 ) ) . < DET I ( 1 ) , DUM (206) ) 

EQUI VALENCE (PIRNEr( 1),DJM(220)),(PIROlD(1), DUM (233) ) 

EQUI VALENCE(PI lNbR( 1) ,DJM(246) ) , (PI lOLD( 1) ,DUM(260 ) ) 

EQUIVALENCE ( AUMR9 { 1 , 1 ) . DUM ( 273 ) ) , ( AUM 19(1,1), DUM ( 386) ) 

EQUI VALENCE ( AUMR1 2 ( 1 , 1 ) , DUM ( 5 0 0 ) ) , ( AUM 1 1 2 ( 1 , 1 ) , DUM ( 655 ) ) 
EQUlVALENCt(AK9(l),3R(l)),(A<l2(l),GR(10)),(SJGMA(l),GR(25>) 

EQUI VALENCE (GR< 35) , XI 9< 1) ) , ( = T 9 ( 1 ) , GR ( 45 ) ) . ( X U2 < 1 ) , GR ( 6 0 ) ) 
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FORMULA FOR GR , G I , A LP H A , BE T A , G AMM A , DEL T IS 3<ND+1)**2 

MAINFRAME PROGRAM FOR THE NASA FLUTTER 

REAL MU 
N U M M A = 1 
INP = 0 
MOs 9 

LK = NF A SPEC I IF I E 0 INTEGER USUALLY LESS THan 3 


N I N = 5 
NP = 6 
N 0 U T s N P 

1001 FQRMAT(1H0*START*) 

WR I Tt ( NP, 1001 ) 

WR I TP ( NP, 10 01 ) 

L I =3 
L J = 3 
MU= , 9 
GAMMAP=2. 

NF s 3 
LK = NF 
N s NF 
K = 1 

LM = 3 
I. J.t s LJ + 1 
7/31 FORMAT (8F10. 4) 

7345 FORMAT ( 4E13 . 4 ) 

7734 FORMAT ( 1 X , I3**,*l3*,+I3, 4E20,8) 

L=3LLG0 

L=3LMKK 

DO 1575 MM=1,250 

DUM(MM) = 0 
1575 CONTINUE 
KEN=3LK6N 
L2(i)=2LSA 
L2(2)=4LCLAD 
l.2(3)=3LQSF 
L 2 ( 4 ) = 0 

CALL SEGMtNTIKEN, 1, ,2, 1, NJMMA ) 

FOR A FLUTTER Run Call SA for THE AR, AI, BR, bi 

CALL SA(AR,AI,8R,BI,lI,LJ,NP,GAMMAP,MU,DEG) 
L2(1)=3LS1B 
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l 


L2<2)=0 

CALL SEGM£NT(KEN»l#i_2*liO) 

GET C«(IJK) AND CKIJK) 3Y THIS CALL To GET INPUT FOR SECT2A 

CALL S18(OUTPUT,CR,0I, I MPdT, AR, AR » A I . A I # BR . B I , L I • L J» LK. NF , M» ) 

GET THE L)R(IJK), A D D I ( 1 JK > BY THIS CALL 
WHFRE DUM IS A OUmMY SECTOR OF LENGTH BR(IJK) BUT ZERO EVERYWHERE 

CALL S1B( CHJT,DR, UI , IN, AR, 3R. A I # 8I , DUM, IJUH.LI #LJ, LK.NF,NP) 
L2(1)=3LS1C 
L2(2)=D 

CALL SEGMtNT ( KEN, 1 , _2, 1 , 0 ) 

GET THE £H ( I JK ) , AND EI(IJK) BY THIS CALL 

CALL S1CTAR, A I ,CK,C I , CR,CI , DR, UI.ER, El ,LI,LJ,LK,NF) 

GET THE KR(IJK), AND F I ( I jk ) BY THIS CALL 

CALL S1C(BR,B1,CR,CI,DR,DI,DJH,DUM,FR,FI,LI,LJ,LK,NF) 
p I =3 , 141B926 
DELTA= .0U01/PI 
L2 C 1 )=2LSi 
L 2 ( 2 ) =4LPRGD 
L2(3>=0 
KEN=3LLG0 

C GET THE EVALUATION =>ONTS FOR THE 9TH AND 12TH ORDERS 

CALL SEGMENT(KEN,1,.2,1,NUMMA) 

CALL S3(1,PSI9,ETA9,YR9,YI9,<9»IN,DELTA,CR,DR,NF.9,NP) 

CALL S3(1,PSI1?.ETA12,YR12,YU2,K12,IN,DELTA,ER,FR,NF,12,NP) 

OU T = 0 

FPS= . OOOOUQOOOQOl 
DELTAX= .0000000001 
IN: 0 
L = 3 L H K K 
L 2 ( 2 ) = 2LOP 
L2(3)=3LPKE 

L2(1)=3LS2A 

L2(4)=0 


lit 
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KEN=3LKEN 

CALU SEGMENT ( KEN , 1 , 1.2 » 1 , NUM'IA ) 

D DOES SETUP AND FINOS ALL CONVERGED DETERMINANTS 

CALL D(FFR» FF I ,GR,GI .ALPHA, BETA, GAMMA, DEIT.CR# Cl # DR. DI.AUMR9.AUMI9. 
1,AUMR12,AUMI12 ,dYR9,DYI9»PIR3LD,PIRNEW»YR9.YI9.NDR9,NDI9.9.K9.1,N> 
CALL D(FFR.FFI,GR.6I , ALPHA,9ETA,GAMMA,DELT,ER,EI *FR,FI , AUMR9. AUMI9 
1.AUNR12, AUMI12.DYR12, D Y 1 1 2 , P I I OLD , P I I NEW , YR12 , Y 1 12 , NDR1 2 , ND 1 12 , 
212.K12.1.N) 

L2(1)=2LS4 
L*3LMKK 
L2 ( 2 ) = 0 

SET UP THE 9 TH ORDER SIMULTANEOUS EQUATIONS 
TAKING ADVANTAGE OF COMPLEX -DETERMINANT PAIRS 

CALL SEGMENT(KEN#1»l2»1,NUMMA) 

CALL S4(CC.BfGA, A9.39.C9,XR,XI ,U»PSI9,ETA9.9.K9,PI .DYR9.DYI9. YR9. 

1 Y I 9 ) 

SOLVE THE SIMULTANEOUS EQUATIONS 
CALL SIMQ(CC,9,BIGA) 

SET UP THE 12TH ORDER SIMULTANEOUS EQUATIONS 

CALL S4(CCH,BIGAH, AH, 3H, CH, XRH,X I H.U. PI S12.ETA1 2* 12.K12.PI » 

1 DYR12,DYH2, YR12.YU2) 

L2(1)=4LS1MQ 
L2 ( 2 ) = 0 

CALL SEGMENT<KEN_.1»l2.1.NUMMA) 

SOLVE THE SIMULTANEOUS EQUATIONS 

CALL SIMQ(CCH,12#3I3AH> 

1106 L2 ( 2 ) = 0 

L2(1)=3LS5A 
L 2 ( 2 ) = 3LTE A 
L2(3)=4LPRQD 
L 2 ( 4 ) = 0 

CALL SEGMENT < KEN< 1 » L2 * 1 , NUMM A ) 
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SET UP THE POLYNOMIAL COEFFICIENTS For THE 9TH ORDER AND GET 

CALL S5A(0UT,AK9» IN,K9,4,BIGA,B9,C9, IERR) 

1177 A Y ( 1 ) = AK9 ( 9 ) 

A Y ( 2 ) s AK9 ( 8 ) 

A Y ( 3 ) = AK9 ( 7 ) 

AYM) * AK9(6) 

A Y ( 5 ) s AK9 ( 5 ) 

A Y ( 6 ) = AK9(4) 

A Y ( 7 ) s AK9 ( 3 ) 

A Y ( Q ) * AK9 ( 2 ) 

A Y ( 9 ) = AK9( 1 ) 

A Y < 1 0 ) = 1, 

GET THE POLYNOMIAL ROOTS 
CALL PRQD(AY,9,RR,RRI, YY,NUM, I ERR) 

THE psi and ETA From SU3R0UTINE tea 


1175 CALL TEA(9,6,RR,RRI,XI9,ET9) 

K12 = KTWELV 

C SET UP THE 12TH ORDER POLYNOMIAL COEFFICIENTS 
C 

CALL S5A(OUT, A«l2, I N, <12, 6 , 3 I G AH , 3H , AH , CH * I ERR ) 
1171 A Y ( 1 ) = AK12 ( 12 ) 

A Y ( 2 ) = A K 1 2 ( 1 1 ) 

A Y ( 3 ) = A K 1 2 ( 10 ) 

A Y ( 4 ) = A K 1 2 { 9 ) 

A Y ( 5 ) = AK1 2 ( 8 ) 

A Y ( 6 ) = AK12(7) 

A Y ( 7 ) = AK12 ( 6 ) 

A Y ( 8 ) s AK12 ( 5 ) 
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A Y ( 9 ) s AK12 ( 4 ) 

AY(10) = AK 1 2 ( 3 ) 

A Y ( 11 ) = AK 1 2 ( 2 ) 

A Y ( 12 ) = AK12 ( 1 ) 

A Y ( 1 3 ) = 1, 

GET THE HOOTS FROM PROD 

CALL PRwD(AY,12»R12TR, R12TI.YY,NllM,lERR) 

GET THE 12TH ORDER P SI AND ETA FROM TEA 

1169 CALL TEA(12,6,R12TR,R12TI ,XI12,ET12) 

ON TO SECT 6 
L2(1)=2LS6 
L 2 ( 2 ) = 0 
L2<3) = 0 

KENs 3 L L G 0 

CALL SEGMtNT(KEN*L2.1f NJMi*) 

GET THE COMMON POLYNOMIAL COEFFICIENTS 

CALL S6(0UT,SIGMA,I\|,AK9,AK12,NP) 

1167 Z ( 7 ) = 1 

7(6) = SIGMA(l) 

Z ( 5 ) = S I GM A ( 2 ) 

Z ( 4 ) = S I GM A ( 3 ) 

Z ( 3 ) = S I GMA ( 4 ) 

7(2) = S I GM A ( 5 ) 

7(1) = S I GM A ( 6 ) 

KEN=3LLG0 
L2(1)=4LPR0D 
L2(2) =36TEA 
L2(3)=3LPAT 
L2 ( 4 ) =0 

CALL SEGMENT<KEN>1,L2»1,NlIMMA) 

GET THE 6TH ORDER ROOTS 
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CALL PRQD(Z»6,R06TR,R06TI ,YY,NtjM, IERR) 

GET THE 6 T H ORDER =»SI AND ETA FROM TEA 

CALL TEA(6«6«R06TR»3Q6TI,XI6>ETA6) 
ga=gammap 


GET THE FINAL OUTPUT SHEETS 

CALL PAT<AUMR9#AUHRl2»PSl9,PSI12#Xl9»XU2,Xl6,eTA9.ETAl2.ET9»ETl? 
1.ETA6, GA,MU, NF» RR» RRI » R12TR.R12TI ,R06TR,R06T I .NDR9. NDR12, YR9, YR 
212»PIR0Ln.PlRN6w»PII0LD,Pl!N = W,AUMl9.AUMll2,NDl9*ND112»YI9,YU2) 
STOP 
END 

NOLIST 

•END 
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SUBROUTINE S A ( AR , A I , BR , 3 I , L I , L J , NF , GAMMAP , MU , DEG > 

DIMENSION Y ( 30 ) * Y Y ( 30>,W< 30>,WW< 30),ZZ< 30>,Zl( 3O),AMO(30), 
1 CLA(30)iAPHU(30),3PMU(30) » AR < 125 > , A I < 125 > , BR < 125 ) , B I < 125 > 
DIMENSION APHI2(30),8PNI2(30) 

SUBROUTINE FOR GENERATING THE AR, Al,8R,8l, FO THE NASA FLUTTER 
EQUATIONS. 

REAL MU,MU1,MUSIPH, MUCOPH, MZZERO , MUZERO, M 0 XPH I , Mo ZPH I 
C SET UP EXTERNALS CLA(I), AMO(I), gammap, mu 
LK *NF 
NF2 * 2 . *NF 
ANF2 = NF2 
DEG* 360./ANF2 
NP = 6 
LJ1 = L J+l 
AMO ( 1 ) = 0 
AMO ( 2 ) = .3 
A MO ( 3 ) = .4 
AMO ( 4 ) = .5 
AMO ( 5 ) = .6 
A MO ( 6 ) = .65 
AMO ( 7 ) = .7 
AMO ( 8 > = .75 
AMO ( 9 ) =1. 

CL A ( 1 ) = 6,398 
CL A ( 2 ) = 6,398 
CL A ( 3 ) =6.517 
CL A ( 4 ) = 6.876 
CL A ( 5 ) = 7,792 
CL A ( 6 ) = 8,594 
CL A ( 7 ) = 9,072 
CL A ( 8 ) = 10.027 
CL A ( 9 ) = 10.027 
DO 50 I “ 1 , L I 
DO 50 J = 1 , L J 
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DO 50 K s 1 i NF 
IJK * ( I*Ll-LJl+J)*LK*K 

AR(IUK) s o 
A I < I UK ) ? 0 
RR(IJK) = 0 
8 1 ( I UK ) * 0 
50 CONTINUE 
NINT s 20 
ANINT * NINT 
G = ( .97". 2)/ ANINT 
H e (l,-,2)/ANlNT 
PI * 3,1415926536 
P I 2 = PI/2, 

NF2 = NF* 2 

ANF » NF 

RNF = l./ANF 

PHI =(PI*ANU)/ANF 

MU1* 1, ♦ MU 

P 1 180 = PI/180, 

NDEG = 360, /DbG 

NINT=NINT+1 

DO 1 NPHI s l.NDEG 

ANPHI = NPHI-1 

PHEYE- ANPHI*DEG 

PHI s ANPHI *P I 180*DEG 

COSPHI = COS (PHI) 

SINPHI = S I N ( PH I ) 

MUSIPH = MU*SINPHI 
MUCOPH s MU*COSPHJ 
NINT1 = NINT+1 



no 2 J= 1 » N I N Ti 
CMAY = J-l 
L = J + l 
Z = ,2+CHAY*G 
X s ,2+ C H A Y * H 
2 MU = I* MUSIPH 
EMU = X + HUS 1 PH 
HZZERO= ABS(ZMU) 

MUZEHO = ABS(fcMU) 

IF(EMU) 52.53,53 
53 AQXPHI = 0 
GO TO 55 
5? AQXPHI = pi 
55 MQXPH I = ( ,8*MUZtR0) /HU1 
IF(ZMU) 56,57,57 
57 AOZPHI = 0 
GO TO 59 

56 AOZPHI = PI 
59 MOZPHI = ( .8*MZZER0)/MU1 

CALL CLAD(PI,AMO,CLA,AOXPHI,HOXPHI.CLX,CDX,NP,IER) 
CALL CLAD(PI , AMO.CLA, AOZPHI , HOZPHI .CLZ.CDZ.NP, IEPR) 
Y Y ( J ) = MUZERO*(X-,C5)*(X-.05)*CDX 
Y ( J ) = MUZERO*(X-,05)*COX 

WW(J) = MZZERO*Z*Z*CLZ 
W(J) - MZZERO*Z*ClZ 
ZZ(J) = MUZERO*X*X*CDX 

Z1(J) s muzero*x*cdx 
? CONTINUE 
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DO 2 J= 1,NINT1 
CHAY = J-l 
l = J + l 
Z * , 2 + CHA Y*G 
X s i 2* CHAY*H 
ZMU = Z* MUSIPH 
EMU = X + MUSIPH 
MZZERO* ABS(ZMU) 

MUZEPO = ABS(bMU) 

IF(EMU) 52,53,33 
53 AOXPHI = 0 
GO TO 55 
5? AQXPHI = PI 
55 MQXPHI = ( ,8*MUZtR0)/MUl 
IF(ZMU) 56,57,57 
57 AOZPMI = 0 
GO TO 59 

56 AOZPHI = PI 
59 M0ZPH1 = ( , 8*MZZbR0 ) /MU1 

CAUL CLAD(PI,AMO*CLA, AQXPHI, MOXPHI,CLX, COX, NP.lEW) 
CALL CLAO(PI,AMO,CLA,AQZPHI,'10ZPHI#CLZ,CDZ,NP,IERR> 

YY(J) = MUZERO* ( X- , Q5 ) *( X- . 05) *CDX 
Y ( J ) = MUZERO* ( x- , 05) *CDX 
WW(J) = MZZ£RO*Z*Z*CLZ 

W(J) = mzzero*z*clz 
Z Z ( J ) = MUZERO*X*X*CDX 
Z1 ( J ) = MUZtRO*X*COX 
2 CONTINUE 


CALI, QSF(H, YY,YY»NIMT) 

CALI QSF(H> Y, Y,NlNT> 

CALL QSF ( G,WW,WW,NINT> 

CALL QSF(G,ZZ,ZZ#NINT> 

CALL QSF (G»W»W*NlNT) 

CALL QSF(H#Zl,Zl,MlNT) 

APHJl(NPHI) = YY(NINT) 

BPH 1 1 ( NPH I ) = MUCOPH* Y ( N I NT ) 

A PH 1 2 ( NPH 1 ) s 0.4286874*(HW(NINT) «• ZZ(NlND) 
8PH 1 2 ( NPH I ) = 0 . 4286874* ( W(VINT) + Zl(NlNT)) 
1 CONTINUE 

DO 15 K * 1 * NF 

Kll * (LI-LJ1+1)*LK+K 

K22= (LI*2-LJ1+2)*L<+K 

A K 1 = K - 1 

PIAK.I s PI*AKl 

VALUE1 = 0 

VALUE2 = 0 

VALUE3 = 0 

VALUE3 = 0 

VALUE4 s 0 

VALUE5 = 0 

VALUE6 = 0 

VALUE7 = 0 

VALUE8 = 0 

HO 10 L : 1 » NF2 

ANU * L-l 

PIK1NU = (PI AK1*ANU)/ANF 
COSPIsCOS(PlKlNg) 

SINPI = SIN(PIKlNU) 

VALUE1 s VAlUEI +AP4I1(L)*C0SPI 
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VALUE3 s VALUE3 ♦ BPRIKU • COSPl 
VALUE4 s VALUE4 *BP4 1 1 < L ) *S I NP I 
VALUE5 s VALUE5-t-ApHI2(l)*C0S»I 
VALUE6 s VA|_UE6 + APHI2(L)*SIN 9 l 
VALUE7 = VALUE7 ♦BP4I2<L> *COSPj 
VALUES s VALUE8 + B»HI2(L>* S X Np I 
10 CONTINUE 

AR(Kll) = VALUEl*RNP*3A^MAP 
AR ( K22 ) * VALUE5*RNF *GAMMAP 
AI(Kll) s -VALU62*RNP *GAMMAP 
A I ( K22 ) * -VALUE6*RNP*GAMv|AP 
RR(Kll) = 2 , *VALUE3*RNF ♦ GA'IMAP 
HR ( K22 ) * 2. *VALUE7*RNF ♦GAMMAP 
HI(Kll) = -2,*VALUE4*RNF *GAMMAP 

R I ( K22 ) * -2,*VALUE3*RNF *GAMMAP 

15 CONTINUE 

K221 =(2*LI-LJ1+2)*lK+ 1 
BR ( K221 ) = BR(K22l>*4. 

Kill =(LI-LJ1 + 1)*LK+1 
BR ( Kill ) = BR ( K 1 1 1 ) * , 3157896 
K133 = (3*LI-LJ1+3)*LK+1 

K233 = (3*LI-LJl+3)*L<+2 

AR(K133)=2. 

BR ( K133 ) = 2, 

ARCK233) = 1,0 

RETURN 

END 
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SUBROUTINE CLAD(Pl,AM.CL,A,AMO.CLX,CDX,NP,IERR> 
DIMENSION A M ( 1 ) , Ctd) 

I ERR = 0 

IF < A-Pl >1,2, 1 

2 CLX* 5.73 
CDX= 0.02 
RETURN 

1 I F ( A ) 999.4,999 

4 C D X = 0 , 0 1 

3 MA s A MO*l 0 . + 1, 

IF(MA-ll) 40,40.41 

41 IERR=2 
RETURN 

40 GO TO ( 100, 101, 101,101 , 104, 105. 106. 107, 108, 108. 108 ). MA 

100 CLA = CL ( 1 ) 

GO TO 111 

101 CLAs ( AMO- AM ( 2)>*< Cl < 2 > -CL < 1 > ) +CL < 2 ) 

GO TO 111 

104 OLA = ( AM 0 - AM ( 3 ) ) * ( CL ( 4 ) -CL < 3 ) ) + CL<3> 

GO TO 111 

105 CLA = ( AM0-AM<4) )*(CL(5)-CL(4) )+ CL(4) 

GO TO 111 

106 IF(AM0-AM( 6 ) ) 6 0 , 6 0 , 61 

60 CLA s (AM0-AM(5)> * ( CL < 6 ) -CL ( 5 ) ) ♦ CL<5> 

GO TO 111 

61 CLA s (AMQ-AM(6)>* < CL ( 7) -CL< 6> >♦ CL<6) 

GO TO 111 

107 IF(AM0-AM( 8) ) 63* 103,103 

63 CLA =(AMO-AM(7)) * ( Ct ( 9 ) -CL ( 7 ) ) + CL<7) 

GO TO 111 

108 CLA = CL ( 8 ) 

111 CLX=CLA 

RETURN 

999 I N a 1 
RETURN 
END 
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SUBROUTINE QSF < H» Y * Z » ND I M ) 

QSF 

041 


DIMENSION Y(1),Z<1> 

QSF 

044 


HT»,3333333*H 

QSF 

046 


INTEGRATION BY SIMPSONS RULE 
1F<N0IM-5)7,8,1 

QSF 

047 

1 

SUM1sY(2)+Y(2) 

QSF 

050 


SUMlsSUMl+SUMl 

QSF 

051 


SUM1*HT*( Y(l)+SUMl*Y(3> ) 

QSF 

052 


AUX1=Y(4)+Y(4) 

QSF 

053 


AUX1*AUX1+AUX1 

QSF 

054 


AUX1=SUM1+HT*( Y(3)+AUX1*Y(5) ) 

QSF 

055 


AUX2=HT*<Y(1)+3,875*(YC2>+Y(5) ) +2 . 625 * ( Y ( 3 >* Y ( 4 ) )+Y(6) ) 

QSF 

056 


SUM2=Y(5)+Y(5) 

QSF 

057 


SUM2*SUM2+SUM2 

QSF 

058 


SUM2*AUX2-HT*(Y(4)+5UM2+Y(6) ) 

QSF 

059 


7<1)=0. 

QSF 

060 


AUXsY ( 3 ) * Y ( 3 ) 

QSF 

061 


AUXsAUX+AUX 

QSF 

062 


Z(?)=SUM?~HT*< Y(?)*AUX+YC4> ) 

QSF 

063 


7(3) = SUM1 

QSF 

064 


Z(4)=SUM2 

QSF 

065 


IF(NDIM-6)5,5,2 

QSF 

066 



QSF 

067 


INTEGRATION LOOP 

QSF 

068 

? 

DO 4 I = 7 * N D I M , 2 

QSF 

069 


SUM1= AUX1 

QSF 

070 


SUM2= AUX2 

QSF 

0 71 


A UX 1 s Y ( I - 1 ) ♦ Y ( 1-1) 

QSF 

072 


AUX1=AUX1+AUX1 

QSF 

0 73 


AUXl=SUMl+HT*(Y(I-2)+AUXl+Y(I)) 

QSF 

074 


Z( I «? ) =SUM1 

QSF 

075 


I F ( I*NDIM)3,6,6 

QSF 

076 

3 

AUX2 = Y( I ) + Y ( I ) 

QSF 

077 


AUX2* AUX2+AUX2 

QSF 

073 


AUX2=SUM2+HT*< Y< 1-1 ) +AUX2 + Y< I+D > 

QSF 

079 

4 

Zl !-D=SUM2 

QSF 

080 

5 

Z(NOIM-l)=AUXl 

QSF 

081 
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Z(NDIM)=AUX2 

QSF 

08 


RETURN 

QSF 

08 

6 

Z { ND I M-l > = SUM2 

QSF 

08 


Z ( ND I M ) = AUX1 

QSF 

08 


RETURN 

QSF 

08 


END OF INTEGRATION LOOP 

QSF 

08 



QSF 

08 

7 

IF(NDIM-3)12,11,8 

QSF 

08 



QSF 

09 


ND I M IS EQUAL TO 4 DR 5 

QSF 

09 

8 

SUM2sl,l25*HT*(Y(l)*Y(2)+Y<2)+Y<2)+Y<3)*Y<3)+Y<3>+Y(4) ) 

QSF 

09 


SUM1=Y(2)+Y(2) 

QSF 

09 


SUHlsSUMl+SUMl 

QSF 

09 


SUMlsHT*( Y(1)+SUM1+Y(3> ) 

QSF 

09 


Z ( 1 ) s 0 , 

QSF 

09 


AUXlsY(3)+Y(3) 

QSF 

09 


AUX1*AUX1+AUX1 

QSF 

09 


Z(?)=SUM2“HT*(Y(2)+AUX1*Y(4)> 

QSF 

09 


IF(NDIM-5)10*9»9 

QSF 

10 

9 

AUX1*Y(4)+Y(4) 

QSF 

10 


AUX1* AUXl+AUXl 

QSF 

10 


Z(5)=SUM1+HT*( Y(3)+4UX1*Y<5) > 

QSF 

10 

10 

Z(3)=SUH1 

QSF 

10 


Z(4)sSUM2 

QSF 

10 


RETURN 

QSF 

10 



QSF 

10 


ND I M IS EQUAL TO 3 

QSF 

10 

11 

SUM18HT*(1,25*YC1)*Y(2)+Y(2)-,25*Y(3) ) 

QSF 

10 


SUM2=Y(2>+Y(2) 

QSF 

11 


SUM2=SUM2+SUM2 

QSF 

11 


Z(3)sHT*(Y(1)+SUM2 + Y(3>) 

QSF 

11 


Z ( 1 ) s 0 , 

QSF 

11 


Z(2)=SUH1 

QSF 

11 

12 

RETURN 

QSF 

11 


END 

QSF 

11 
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SUBROUTINE SIB ( OUT , CR » C l , I N . AR# BAR » A I , BA I , BR , 8 1 * L 1 , LJ. LK . NF , NP ) 
DIMENSION CII1), CR(1>, AR(1 ), BAR ( 1 ) , BAUD, Al(l), BR(1>, 81(1) 


CALL SECT1(C«,CI , AR, AR. AI, AI , 3R # B I , L I » L J , LK * NF ) 

FOR THE CALCULATION OF THE CR<I.J»K) AND THE CI(I,J,K) USE.., 

FOR THE CALCULATION Or THE DR(I,J,K> AND THE Dl(I,J,K) USE.., 

CALL SECT1(DR,DI , AR,3R, AI ,BI ,DUM,0UH,LI .LJ*LK#NE) 

WHERE DUM IS A VECTOR AS LONG AS BI AND BR BUT = TO ZERO EVERYWHERE 

THE IMPUT TO THIS PROGRAM CONSISTS OF THE COMPLEX FOURIER COEFFICIENTS 
OF 18 PERIODIC FUNCTIONS. SPECIFICALLY# THE PROGRAM IS TO Bg SUPPLIED 
WITH. , . 

AREAL ( I # J# K ) AJMAG(I,J,K) 

BREALU»J»K> BIMaG(I,J,K) 

FOR I = 1# 2, 3 
FOR J = 1, 2, 3 
FOR K = 1 » 2 » . , . # NF 

WHERE NF IS A SPECIFIED INTEGER 

THIS SUBROUTINE IS A BRANCH THAT PROCEEDS FROM THE CALCULATION 
OF COMPLEX FOURIER COEFFICIENTS FOP 18 FUNCTIONS RELATED TO THE ABOVE, 


EQUATIONS FOR COMPUTATIONS 


LJ = LJ 

+ 1 


LJ1 = 

LJ- 

i 

LK = NF 



N2F1 

= 2 

* NF - 

NF 1 * 

NF 

♦ 1 

KL * 

2*NF 

-1 

DO 5 

1=1# 

LI 

DO 5 

J = 

1# LJ1 
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c 


IKL = <L1*I-LJ + J>*Kl 

IJ=(U * I «LJ+J)*LK 
C FOR K = 1» 2* ... > Ni- 

C 

DO 55 K = 1 , MF 

AK = K-l 
IJKL = IKL + K 
UK = IJ + K 

CI(IJKL) = 2.* AK * BAR(IJK) ♦ B l C 1J« ) 

CR(IJKL)= -2, *AK*3AI ( UK) + 3H<IJK> 

C PERFORM THE SUMMATIONS 
CONSTR =0 
CQNSTI = U 
C 

DO 7 1= 1, 3 
II = (LI*I- LJ ♦ L ) *LK 
UL S < L I * L -LJU)*L< 

CONST1 =0 
CONST 4 =0. 

C 

DO 15 N = 1 , K 

ILN * IL*N 

LJK1N = JL «• K +1 -N 
ARILN = A K C I L,N ) 

RAIJK1=BAR (LJK1N ) 

A I I LN = A I (UN) 

BlLJKlsBAI (LJK1N) 

APPLE= BALJK1*ARILN 
PEAR = Al ILN*BILJK1 
CONST1 = CONSTl + APPi.E-P6AR 
C C0NSTlsC0NSTl^ARlLN*8ALJKl-Al ILN*BlLJKl 

apple=ariln*biljki 

C CONST 1 = CONST1 * A R ( ILN) * 3AR(LJKlN) - Al(lLN) * BAI(LJKIM) 

C C0NST4 s CONST 4 + A R ( I L N ) * SAl(LJKlN) + A I < 1 L N ) * 8AR ( L JK1N ) 

PEAR s AI1LN+BALJK1 
C0NST4=C0NST4+ APPLE +PEAR 
15 CONTINUE 
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noon 


I 


c 

C NOTE THE LAST SUM ON N IS OMITTED FOR K = NF 

C0NST2 = 0 
CGNST5 = 0 
IF(NF-K) 17, 17, 16 
16 NFK*NF-K 
C 

DO 20 N= 1, NFK 
JLN1 = I L * N ♦ 1 
t JNK= JL+N+K 
ILNK = IL ♦ N + K 
LJN1 = JL ♦ N +1 
ARILN1 = AR( ILN1) 

BRLJNK=BAR<LJNK) 

AIILNlsAl ( ILN1) 

BILJNKsBAl (LUNK) 

ARILNK=A«( ILNK) 

BRLJN1=BAR<LJN1) 

AI I LNK = A I ( ILNK) 

BILJN1=BAI (LJN1) 

APPLE=ARILN1*6RLJNK 
PEAR = AI ILN1*BILJNK 

peach=arilnk*brljni 

PLUM = AI ILNK*BILJN1 

CQNST2 = CQNST2+ARPLE+PEAR+PEACH+PLUM 
1 BILJN1 

C0NST2 = C0NST2+ARILN1*BRLUNK + A I ILN1*BILJNK + ARILNK*BRLJN1^A I IL-NK* 
C0NST2 = C0NST2 ♦ AR < l LN1 ) * BARUJN*) ♦ AldLND * BAKLJNK) + 
1 AR(ILNK) * BAR(LJNl) + AlULNK) * B A I ( L JNl ) 
APPLE=ARILN1*BILJNK 
PEAR = A I ILN1*BRLUNK 
PEACH = ARlLNK*BlLJNl 


63 



non n n n on 


PLUM s A 1 I LNK*BRL JN1 
C0NST5 = C0NST5 ♦ APPLE-PEA R-PEACH*PlUM 
CQNST5 = CQNST5*APILN1*8ILJN<-AIILN1*8RLJNK-ARILNK*9!LJNJ>AI ILNK* 
1 BRLJN1 

C0NST5 = C0NST5 * AR< ILM1) * BAKLJNK) - A I < I LNl > * 8AR(LJN<) - 
1 AR(ILNK) * BAMLJMl) * A I ( ILNK ) * BAR(LJNl) 

20 CONTINUE 

17 CONSTR= CUNSTR «• CONSTI + C0MST2 
CONSTI = CONSTI ♦ CONST 4 + C0NST5 
7 CONTINUE 

CR ( I JKL ) a CR(IJKL)- CONSTR 
CI(IJKL)* CKIJKD- CONSTI 
55 CONTINUE 

FOR K = NF+1, nF + 2, ,,,, 2NF-1 

DO 9 K = NF 1 , N2F 1 
CONSTRsO , 0 
CONSTI = 0,0 
I JKL = I KL + K 

e 

DO 8 L= 1,3 
IL= < L I * I -LU+L)*LK 
JL = < L I *L -LJ+J)*UK 
CONST3 = 0.0 
CONST 6 = 0,0 
KIN s K * 1 - NF 
C 

2? no 30 n = kin ,nf 
ILN = I L+N 

LJK1N = JL + K + 1 -N 
ARILN=AR( ILN) 
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BRLsJK1 = BAR<LJK1N> 

AI I L N = A I ( ILN) 

BILJK1=BAI (UJK1N) 

C0NST3 - C0NST3+ARIuN*BRLJK 1-AI IUN*BlLJKl 
C0NST6+ARILN*BIIJK1*AIIIN*BRLJK1 


C0NST6 = 
C0NST3 = 
C0NST6 = 
CONTINUE 


CONS T 3 ♦ AR(IUN) * BAR(LJKlN) - Al( 
C0NST6 * AR(ILN) * 9AKLJK1N) ♦ Al( 


ILN) 

ILN) 


BAI (LJK1N) 
BAR(LJKIN) 


CONSTR = CONS TR ♦ C0NST3 
CONST I =CONST I + CONST 6 
CM 1 JKL)»-CONSTl 
CR(IJKL)= - CONSTR 
CONTINUE 


CONTINUE 

LJ= LJ-1 
RETURN 
FEND 
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SUBROUTINE PROD < XCOF , M , ROOT R, ROOT I , COF , NUMBER) 

DIMENSION XCOF < 1 ), C OF < 1 ), ROOTR <1), ROOT I (1) 

DOUBLE PRECISION X Q , XO » Y 0 , YO » X * Y , XPR » YPR » UX , U Y , V » YT , XT . U 
DOUBLE PRECISION XT2, YT2, SUMSQ, DX. DY» TEMP, ALPHA 
COMPUTES THE REAL AND COMPlEX ROOTS OF A POLYNOMIAL 
USING THE NEWTON RAPHSON ITERATION TECHNIQUE 
C PARAMETERS 

C XCOF VECTOR OF M + l COEFFICIENTS OF THE POLYNOMIAL 

c ordered from smallest to largest power 

C COF A WORKING VECTOR OF SIZE M + l 

C M THE ORDER 0" THE POLYNOMIAL 

C ROOT I RESULTANT VECTOR OF LENGTH M OF IMAGINARY PARTS 

C ROOT I ( 1 ) IS THE INITIAL VALUE OF THE Y GUESS IMAGINARY PART 

C ROOTR RESULTANT VECTOR OF LENGTH M OF REAL ROOTS 

C ROOTR(l) IS THE INITIAL VALUE OF THE Y GUESS REaL PART 

C IER ERROR CODE 

C IER = 0 NO ERROR 

C SUBROUTINE POLR T ( XC OF , COF , M , ROOTR » ROOT I , I ER ) 

C IER =1 M IS LESS THAN 1 

C IER = 2 M IS GREATER THAN 36 

C IER = 3 UNABlE TO DETERMINE ROOTS IN 50 ITERATIONS 

C IER = 4 HIGH ORDER COEFFICIENT IS ZERO 

C PROGRAMMED BY KG 8LEMEL 

C R , A , S » A , INC. 

C ROCHESTER N.Y, 716 I 

C ROCHESTER N , Y. 716 27l 3450 

C SUBROUTINE POLRT ( XCOF , COF , M , ROOTR . ROOT 1 , l ER ) 

IFIT □ 0 
N = M 
IER =0 

IF (XCOF ( N + l ) ) 1 0 , 25 * 10 
10 l F ( N ) 15 # 15 » 32 
15 IERsl 
20 RETURN 
25 IER=4 
GO TO 20 
30 IER=2 
GO TO 20 
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32 I F ( N'*36 ) 35,35,30 
35 NXsN 

NXX * N + l 
N2 = 1 
KJ1 = N + l 
DO 40 L*1,KJ1 
MT s KJ1«L>1 
40 COF(MT) * XCOF ( L > 
45 XO * .00500101 
YO s 0,01000101 
IN = 0 
50 X a XO 

XO = -10,*YO 
YO = -10, *X 
X s X 0 

V = YO 

IN * IN *1 
GO TO 59 
55 IF I T =1 
X PR s X 
YPR = Y 

59 ICT = 0 

60 UX = 0,0 
UY = 0,0 

V a 0 . 0 
YT = 0,0 
XT a 1.0 

U a COF ( N + l ) 

I F ( U ) 65,130,65 
65 DO 70 I s 1 , N 
L a N - I +1 
TEMP a COF ( L, ) 

XT2 = X*XT-Y*YT 
YT2 * X * Y T +Y*XT 

V a V + TEMP *YT2 
U a U + TtMP*XT2 
F I a I 

UXa!JX-*-F I +XT+TEMP 
UY=UY-FI*YT*TEMP 
XT = XT2 
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70 YT = YT2 

SUMSQ s UX*UX+UY*UY 
JF(SUMSQ) 75,110# 75 
75 DX»( V*UY-U*UX )/SUMS3 
X * X ♦ DX 

DY s - <U*UY «• V*UX)/Sl) s 1SQ 

Y s Y + DY 

78 IF(DA8S(DY)+DABS(DX)-1,D-12) 100,80,80 
80 ICT * ICT +1 

IFUCT-500 ) 60,85,8: 

85 IF(IFIT) 100,90,100 
90 IF < IN — 5) 50,95,95 
95 IER = 3 
RETURN 

100 DO 105 L = 1, NXX 

MT = KJ1 -L*l 

TEMP s XCOF(MT) 

XCOF(MT) = COP < L > 

105 COF(L) = TEMP 
I TEMP = N 
Ns NX 
NX = ITEMP 
IFCIFIT) 120 , 55,120 
110 IF(IFIT) 115,50,115 
115 X * XPR 

Y s YPR 
120 IF IT = 0 

I F ( X ) 122,125,122 

122 IF(OABS(Y/X)-1,D-10) 135,125,125 
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125 ALPHA * X + X 

SUMSQ s X*X + Y*Y 
N s N-2 
GO TO 140 
130 X = 0,0 

MX s NX -1 
NXX 5 NXX -1 
135 Y s 0 . 0 

SUMSQ = 0.0 
ALPHA = X 
N s N - 1 

140 COF ( 2 ) = COFC2) ♦ ALPHA* COFU) 

KsL 

145 DO 150 L = 2 * N 
K s K 

150 COF ( L + l ) = COF ( L*1 > *ALPHA*COP (L ) - SUN|SQ*COF ( L*1 > 
155 ROOT I ( N2 ) = Y 
ROOTR ( N2 ) = X 
N2sN2+l 

IF(SUMSQ>160, 165*160 


160 

Ys-Y 

SUMSQ 

= 0 . 


GO TO 

155 

165 

1 F ( N > 

20*20,45 


END 
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SUBROUTINE SIC ( A R , A I , C R , C I , ACR , AC I , Dr , D I , ER , E I » L I , LU , LK , NF ) 
DIMENSION CR(1),CI(1),DR(1),DI<1),ER(1),EI<1) 

DIMENSION AR(1), Al(l), A CR ( 1 ) , A C I < 1 ) 

ROCHESTER APPLIED SCIENCE ASSOCIATES 
100 ALLENS CREEK ROAD 

ROCHESTER, N, Y, 716 271-3450 

THIS SUBROUTINE IS FOR THE C0M B UTATI0N OF ER(I,J,K) AND EI(I,J,K) 
TO COMPUTE ER(I,j,k) AND EI(I,J,K) USE,,. 

CALL SECT1CIAR, Aj ,CR,CI ,CR,CI ,DR,DI ,ER,EI,LI,LJ,LK,NF) 

OR THE CALCULATION OF FR<I,J,K) AND Fl(I,J,K) 

CALL SECT1C(8R,BI,CR,CI,DR,DI,DUM,0UM,FR,FI,LI,LU,LK,NF) 

LK IS THE LENGTH OF < 

LK=3*NF-2 
LK?NF = 2* NF - 1 
LKNF = NF 
LJ=LJ+1 
LJ3=LJ+1 
LJ1 = LJ-1 
DO 5 1=1, LI 
DO 5 J= 1,IJ1 
IJ = (LI*I-LJ*J)*LK 
IJ2 = (LI*I'LJ+J)*L<2NF 
DO 45 K = 1 , NF 
AK= K + K-2 
I JK= I J+K 
I JK2 = IJ2 + K 

ER(IUK) = -AK-ACI ( I JK2)*DR( I JK2) 

EHIJK) = AK*ACR ( IJ<2)+DI ( I J<2 ) 

C ERC I JK)=-AK*ACI ( I JK)+DR( I JK ) 

C E 1 ( l JK ) = AK*ACR( I J.O-fDI ( I JK) 

CONSTR=0 
CONST I = 0 
DO 55 L = 1 , 3 

IL = (LI*I-LJ+L)*LK2NF 
JL = (LI*L-LJ-*-J)*LKNF 
C JL =(LI*L-LJ+J)*LK 

C IL = ( L I * I -L J + L ) *L< 

CONST 1 = 0 
CONST 7 = 0 

nfk=nf-k 

IF(NFK) 3,3,4 
4 DO 15 N = 1 , NFK 


70 



ILN1*IL+N+1 
L JKN* JL+K+N 

C0NST1=C0NST1+CR< ILV1) *AR<LJ<N> *C I < I LNl ) *A I ( L JKN ) 

C0NST7= CR(IUNl) • Al(LJKN) - C l< I LNl > *AR <L JKN > * CONST 7 
15 CONTINUE 
3 CONST2=0,0 
C0NST8*0 
NFlsNF-1 
DO 20 N= 1# NF1 
ILKNb il+k+n 
LJNlsUL+N+1 

C0NST2 = CR( ILKN) * ARdJNl) C I ( I LKN ) A I ( L jNl ) + C0NST2 
C0NST8 = C0NST8 * CRULKN)* A I ( L JNl ) + C I ( I LKN ) *AR ( L JNl ) 

20 CONTINUE 

C0NST3 =0 
C0NST9 = 0 
DO 25 Nsl.K 
IL.K1N = IL + K + l-N 
LJN* JL+N 

C0NST3 = CONS T 3 * CR<IU1N) * AR ( L JN ) - CKILK1N) * AI(LJN) 
C0NST9= C0NST9 + CR(ILKIN) • AI(LJN) + Cl(lLKlN) * AR ( L JN ) 
25 CONTINUE 

CONS TR= CONS TR+ C0NST1+C0NST2*C0NST3 
CONST I sCONST7 + C0NST8 + C0NST9 +CONSTI 
55 CONTINUE 

ER( I JK) = ER( I JK>- CONSTR 
El ( I JK)s fcl ( I JK)- C3NSTI 
45 CONTINUE 

NF1= NF+1 

NF21 = 2*NF ”1 

DO 60 K = NF1.NF21 

I JK* I J+K 

IJK2 = IJ2+K 

AK * K «• K -2 



EIUJK) s AK* ACR( I JK2UDI ( I JK2) 

ER(IJK) s -AK*ACl ( I JK2) ♦ DR ( I JK2 ) 

C EIUJK) = AK*ACR( IJO+Dl ( I JK) 

C ERUJK) = -AK*ACl<I JK)*QR(I J<> 

CONSTRsO 
CONST I = 0 
DO 65 1 = 1 , 3 
NN=2*NF-1-K 
IL * (LI*I-U<-L)*lK2NF 
JU * <LI*L-LJ+J)*LKNF 
C IL=(LI*I-LJ"-U*LK 

C JL ? <LI*L-LJ*J)*L« 

CONST 4 =0.0 
CONST 10 = 0 . 0 
IF ( NN ) 66 * 66 , 6 7 
67 DO 70 N c 1 » NN 

ilkn=il«-k+n 

LJN1 = JL + N*l 

CONST10 = -CR( ILKN)* AKlJNl) * CIULKN) * AR(LJND «. CONSTiO 
C0NST4 = CR(ILKN) * ARUJNl) ♦ CI(IL«N) * A I < LJN1 ) * C0NST4 

70 CONTINUE 
66 CONST5=0.0 
CONST 1 1 = 0 
C 

no 80 N = 1 > NF 
ilkin=il+k*i-n 

L JN= JL+N 

CONST 11 = CONS TU * CR (I LK1N > * A I < L JN > *C I U LK1 N > *AR ( L JN > 
C0NST5=C0NST5>CR( IHIN ) *AR< L JN > »C I < ILK1N) *A I < L JN) 

80 CONTINUE 
C 

CONSTR=CONSTR+CONST4+CONST5 
CONSTI = CONSTI+CONST10* CONSTU 
65 CONTINUE 

ERUJK) = ER('IJK) - CONSTR 
El (UK) s El (I JK) - CONS T I 
60 CONTINUE 
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c 

NFS = 2*NF 

NF32 = 3 * N F - 2 

DO 100 K*NF2,NF32 

I JK=I J+K 
I JK2 = IJ2 + K 
CONSTR = 0 
CONST I = 0 
C 

DO 105 L s 1 » 3 

C 11* < 1*LI-LJ+L)*LK 

IL * ( I*L1-LJ+L)*LK2NF 
C JL *(LI*L-LJ+J)*LK 

JL = (LI*L-LJ*J)*lK\|F 
CONST 6= 0 
CONST 12 = 0 
NN s K-2*NF ♦ 2 
I F ( NF- NN) 101, 102, 102 
C 

102 DO 110 N = NN , NF 
ILKINs I L+K+l-N 
LJNsJL+N 

C0NST6 = CR(ILKIN) • A R < L JN ) - CMlLKlN) • Al(LjN) 
C0NST12 * CR(IUKlN) * AKLJN) * CI(ILKIN) • AR ( L JN ) 
110 CONTINUE 
C 

101 CONSTR* CONSTR CONST6 

CONSTI = CONSTI ♦ C0NST12 


105 

CONTINUE 
ER < I JK ) = 

- CONSTR 


EI(IJK) = 

- CONST! 

100 

CONTINUE 


c 



5 

CONTINUE 


C 

LJ=LJ-1 

RETURN 

END 


• END 
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♦DECK , S4 

SUBROUT I Mb S4 < CC* 8 1 3A , A , B , C* XR , X I , U, S I , E, NO, K , P I , DYR, D Y I , YR, Y| ) 

D I MfeNS I ON XR < 13 ) i X I { 1 3 ) 

DIMENSION A(1),b(i) i C(1),CC(1)iYR(1)#YI(1)*DYr(i) > DYI( 1),X(1),E(1) 
DIMENSION BlQA(l), U(l) ,31(1 ) 

WR IT t (6,144b) 

1445 format ( * entry into S4*> 

SETTING UP OP THE LINEAR AlGEBRAIC EQUATIONS 9TH AND 12TH OROER SYSTEMS 
PROGRAMMED FOR D«, P, CRIMI 
ROCHESTER APPLIED SCIENCE 
ROCHESTER NEW Yo«K 
AC716 271 3450 

PROGRAM FOR SETTING UP THE SYSTEMS OF SIMULTANEOUS EQUATIONS 
FOR THE 9TH AND 12TR ORDER SYSTEMS 
INPUT PARAMETERS 
A ( 1 ) IS DUMMY STORAQE 
R ( 1 ) IS DUMMY STORAGE 
C(l) IS DUMMY STORAGE 
MAXIMUM NECESSARY LENGTH OF A»B,C IS 12 
DELTY IS THE PARAmTER TO BE USED JF ThE DENOMINATOR BECOMES SMALL 

CC < 1 > IS IME RETURNED SET FO SIMULTANEOUS EQUATION MATRICES 
C(M,N> IS OF ORDER C ( 9 , 9 > OR C(l2,l2) THEREFORE C(81) OR C < 1 4 4 > 

Y ARE THE INPUT VALUES OF THE DETERMINANT 
DY ARE THE DELTA Y 
X(l) IS dummy storage MAX IS X c 1 2 > 

E(l) IS ETA(J) MAX IS E T A ( 1 2 ) 

U ( 1 ) IS DUMMY STORAGE MAX DIMENSION IS UU2) 

A PROGRAM For SETTING UP THE SIMULTANEOUS EQUATIONS FOR 9TH AND 12TH ORQE 
NOl = NO r 1 
NK*K*2 

DO 5 Ns 2»NK,2 
A(N)sPI*Si (N) 

B(N> = pI*E(N) 

5 CONTINUE 
N1 = NK*1 
NK1 s NK"1 
DO 10 N=N1,N0 
C(N>ePI*SI (N) 

10 CONTINUE 

DO 25 MM 8 1» NKl, 2 
XR(MM)s-PI*YR(MM) 

XI <MH)«-PI*YI <mm) 

XXR=XR(MM) 

XXIsXI (MM) 

DO 15 N=1,NK1 ,2 
NN=N+1 
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BNN=B(NN)+XXI 

XIR=B(NN)-XXI 

AAsA(NN) 

XA A=XXR+A A 
EXAA=EXP(XAA) 

EMXAA=EXP(-XAA) 

C0SHY=(EXAA+EMXaA)*,5 

U(N)« ( ( SInONN) / ( CO SHY -COS ( 3NN ) ) ) + ( S I N ( X I B ) / ( CoSHY-COS ( X 1 9 ) ) ) )*, 5 
15 CONTINUE 

D02 0 N = 2 » NK » 2 
AA = A ( N ) 

XXAA=XXR+AA 

BNaB(N) 

XI8N=XX I+BN 
BNX I =BN-Xx I 
FXXAA = EXP ( XXAA ) 

EMXXAA=EXP(-XXAA> 

SI NHY= (EXXAA-EMXXAA )*. 5 
C0SHY=(HXXAA+EMXXAA)*.5 
7022 FORMAI (* NOP = % 15) 

U(N)s( (SINHY/( COSHY-COS (X ION) ) )♦ (S I NHY/ (COSHY-COS ( BNX J ))))*,? 

20 CONTINUE 

DO 21 N = N1 » NO 
XXCN=XXR+C(N) 

EXCNsEXP < XXCN ) 

EMXCN=EXP< -XXCN) 

SINHY=(EXCN-EMXCN)*.5 
C0SHY=(EXCN+EMXCN)*,5 
U ( N ) = SlNHY/(COSHY-COS(XXI ) ) 

21 CONTINUE 

DO 30 N= 1 » NO 
MN=(MM-l)*NO*N 
30 CC(MN) = U(N> 

25 CONTINUE 

DO 125 MM=2,NK,2 
XXR=-PI*YH(MH) 

XXI =-P I * Y 1 (MM) 

DO 130 N=1#NK1,2 
XXRA=XXR+A(N*1) 

EXXRA=EXP(XXRA) 

EMXXRA=EXP(-XXRA) 

SINHY=(EXXRA«bMXXRA)* .5 
C0SHY=(£XXRA*hMXXRA)*.5 
BNX I *8 ( N + l ) • XX 1 
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X I BN«8 ( N*1 ) +XX.I 

U <N)*(S I NHY/ (COSHY-COS (9NXI ))- <S INHY/ ( COSHY-COS ( X IBN > ) ) )*.9 
130 CONTINUE 

DO 135 N*2,NK,2 
XXAN*XXR*A(N) 

BNXIs8(N)-XXJ 
X I BN*B ( N ) +XX I 

COSHY=(EXP<XXAN)+EX?( -XX AN) )*,5 

U(N>»(-SIN(BNXI )/ (COSH Y *C0S( 3NXJ ) ) *S I N ( X I BN )/( COSHY-COS < XI BN) > )* ,5 
135 CONTINUE 

DO 140 N*N1,NQ,1 
XRCN*XXR*C ( N ) 

COSHY=(EXP(XRCN)*EX?(-XRCN) >*.5 
U <N)«S IN (XXI ) /(COSHY- COS (XXI ) ) 

140 CONTINUE 

DO 145 N*1 * NO 
MN*(MM-1)*NQ*N 
145 CC(MN)*U(N> 

125 CONTINUE 

DO 150 MM* N1 > NOl 
Ms(MM-l)*N0 
XXRs -P I * YR ( MM ) 

DO 155 N*l*NKl,2 
XRAN*XXR*A(N*1) 

COSHY=(EXP(XRAN)*EXP(-XRAN) )*.5 

MN*M«-N 

BN1* B(N*1) 

CC(HN)«SIN(8N1>/<C0SHY-C0S(BN1) ) 

155 CONTINUE 

DO 160 N s 2 # N K * 2 
MNsM+N 

XRAN=XXR+A<N) 

EXRAN=EXP(XRAN) 

EMXRAN=EXP(-XRAN) 

SINHY=(EXRAN-EMXRAN)* , 5 
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C0SHY=<EXRAN*£MXRAN)*.5 
160 CC(MN)sSiNHY/(COSHY-COS(B(N> ) ) 
DO 165 N»N1,N0 
HNaM + N 

XRCNa ( XXR+C ( N ) )*,5 
EXRN*EXP(XRCN) 

EMXRN=EXP( -XRCN) 
SINHY=EXRN-EMXRN 
COSHY=EXRN*EMXRN 
CC < MN ) s COSHY/SINHY 
165 CONTINUE 
150 CONTINUE 

DO 40 N * 1 > NK 1 , 2 
MNaNO*(NO-l)*N 
MN1 a MN*1 
CC(MNl) a 1. 

CC(MN) = 0 . 0 
40 CONTINUE 

DO 70 N = NK > NO 
MNaN0*N01+N 
CC(MN)=1, U 
70 CONTINUE 

DO 50 M=1,NK1,2 
BIQA(M)=DYR(M)-1, 

50 BIQA(M*l)s-UYl (M*l) 

DO 60 M=N1,N01 
60 RIQA(M)=DYR(M)-1, 

R I GA ( NO ) * 0 

RETURN 

END 

♦ END 
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SUBROUTINE S2A(0UT# PTB.FFI .GS.GI, ALPH,BET.BET#GaMM.DELTfPIR,PH » 

1 IN#CR.Cl,ORt DI » YR, YI .NfNn. U l *LJ fNF» H»NP,NO* EPS# DELTAX) 

C TO GET THE REAL AND IMAGINARY 3 ART OF A COMPLEX DETERMINANT 
C SUBROUTINE SECT2B < OJT , FFR , FF I , GR , G I , ALPH , BET . G AMM , DEL T , P I R , P 1 1 , 

DIMENSION CR(1>, CI(l)*OR(l)#DKl)*ER(l)#El(t) 

DIMENSION FFR(t). FFI <1) 

DIMENSION FR(1).FI (1),3R(1),GI (1) 

ALF(FFRI «GRI J.FFI 1 .31 I J, A0) =(FFRI*GRIJ +FF I I *G 1 1 J ) + A8 

BETF(FFRIfGIIJ,FFH,GRIJ.AB)«(FFRI*GlIJ-FFlI*GRlJ)*AB 

FRF< YR» YI »R.CRMM1,DRMM)=YR**3-3,*YR*( <2 ,*R*Y I )**2)+YR*CRMMl+DRMM 

F I F ( YR. YI ^R.CRMMl) = < 2 , *R + Y I ) * ( 3 . * YR*YR- ( ( 2*R + Y I > **2 ) +CRMM1 > 

GRF(YR,Yl»S,CRM,ClM,DRM)sYR*CRM-<2.*S + YI )*CIM+DRM 

GIF( YR. YI »S,CRM,CIM,DIM)syr*CIM+( 2,*S+YI )*CRM+DlM 

FHI (YR.YI.R»CRHM1)»(2.*R+YI)*(4,*YR*(YR*YR-( ( 2 . *R+Y I > **2 ) >*CRMM1) 
FHR<YR,Yl,R,CRMMl,DRMMl)=(<YR*YR-(2.*R+YI)**2)**2)+YR*CRMMl- 
l(4,*R*YR+2**YR*Yl )**2 +DRMM1 
N I N = 5 
Npa6 

ORDER=NO 

IF<0RDER-9, >1,2,1 

1 IF(0RDER-12. >4,3,4 

1005 FORMAT ( 4 1 H 1 I N DELY NEITHER A 9TH OR 12TH ORDER CALL) 

3 NORs 2 
LKs3*NF-2 
GO TO 8 

2 N0R=1 

LK = 2*NF - 1 
8 LJ1 =3*(2*ND+1) 

C LJ1 IS USED FOR 2 DIMENSION INDICES 
LM=L I 
LNsLJ 
L J=L J+l 

LK33 = L I *L J*LK 

MDs-ND 

MZ = 0 
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MMls(LI*M-LJ+M)*LK*l 
CRMM1 s CR(MMl) 

DRHM1 = DR(MHl) 

DO 89 MR=MD,ND 
RsMR 

I = 3* ( MR+ND) +M 
NIRM=3*(ND-MR)*M 
23 GO TO <66, 67), NOR 
C DO 9 TH ORDER EQUATIONS 

66 FFII=FIF<YR,YI,R,CRMM1) 
FFRI=FRF(YR,YI,R,CRMM1,URMM1) 

ALPH( I)=1,/(FFRI*FFRI+FFII*FFII ) 

FRN= FRF<YR,-YI ,R,CRMM1,DRMMI) 

FIN«FIF<YR,*YI»R»CR'1Hi) 

BET(NIRM)=i,/(FRN*FRN+FIN*FIN) 

GAHM<NIRM)=FRN 

DELT(NIRM)=FIN 

FFI < I ) =FF i I 
FFR ( I )sFFRI 
GO TO 89 

67 FFRI=FHR(YR, YI,R,CRMM1,DRM'U ) 
FFI I=FHI<YR,YI,R,CR'1M1) 

FFI ( I )=FFI I 
FFR ( I )=FFRI 

ALPH( I )=1./<FFI 1*FFI I*FFRI*FFRI ) 
FRN = FHR( YR, - YI ,R,CRv)M 1, DRMMl) 
FIN=FHI(YR,-YI,R.CRMM1) 
RET(NIRM)sl. 7<FRN*FRN*F IN*FI N) 
GAMM(NIRM)=FRN 
DEL T ( N I RM ) =F I N 
89 CONTINUE 

DO 5 MS=MD#ND 
S s MS 

S2 = S+S 

MSND = ND-MS 

MSND3 = MSNO+MSND+MSNO 

NDMS=ND+MS 

NDMS3 = NDMS+NDMS+NDMS 
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a o no ci naan 


DO 5 MR * MS,ND 
R = MR 

NDMR = ND-MR 
NDMR3 = NDMR*NDMR+NDMR 
MRND* MR + ND 

MRND3 = MRND+MRND+MRND 
MRMS1 = MR-MS+1 
DO 5 Msl.LM 
I a MRND3+M 

! s 3*<MR+ND)+M 

THE I INDEX REQUIRES ONLY MR AND M 
INDICES ARE IN THREE DIMENSIONS 
NIRMp 3* ( -MR + ND ) *M 
N l RM = NDMR3+M 
BApBET(NIRM) 

FRNEGI JsGAMM(NIRM) 

FINEGI JsDELT(NlRM) 
NIRM1=(NIRM-1)*LJ1 
IRMe 3*(MR+ND)+M 
FFR I *FFR ( I ) 

FFII = FFI(I) 

AB*ALPH ( I ) 

IRM s I 

IRM1 = (IRM-l)*LJl 
DO 5 N*1#LN 
MNs(LI*M'-LJ + N)*LK 

MNRS1=MN+MR-MS*1 
MNRSls MN+MRMSl 
NJSN s MSND3+N 
NJSNs 3* ( -MS+ND ) ♦N 
NEGI J=(NIRM-1 )*lJ1 *NJSN 
NEGIJ = NIRM1+NJSN 
JSN = NDMS3+N 
C JSN* 3*<MS+ND)+N 
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IJs IRM1+JSN 

C IJ»( IRM-1)*LJ1 JSN 

IF(MRMSl-LK) 6177,6177,6117 
6117 GR(|J) s 0, 

Gl ( I J) =0 
GR ( NEG l J > * 0 , 

G I { NEG I J ) *0 , 

GO TO 5 

6177 IF(M-N)2111, 617, 2111 
617 IF(MR-MS) 2111,2112,2111 
2112 GRUJ) s 1, 

Gl ( I J) = 0 
GR1NEGI J)*l, 

G I ( NEG I J ) =0 , 

GO TO 5 

2111 CRMNRS* CR(MNRSl) 

CJMNRSs CI(MNRSl) 

DRMNRSs DR(MNRSl) 

DIMNRSs DI(MNRSl) 

GRI J=GRF( YR, YI ,S#CRMNRS,CIMNRS,DRMNRS) 

Gl I J = GIF( YR, YI ,S*CRMNRS,CIMNRS,DIMNRS) 

TEMPA s(FFRI*GRIJ ♦ FF I I *G 1 I J ) * A8 
G I ( I J ) = (FFRI*GI I J-FFI I*GRI J)*AB 
GR(IJ) = TEMPA 

3111 GRNEGI J=GRF(YR,-YI,S,CRMNRS,CIMNRS,DRMNRS) 

GINEGI J=G IF (YR,-YI , S, CRMNRS, CIMMRS.DlMMRS) 

Gl (NEGI J)=-BETF cFRN=GI J,GINE3I J.FINEGI J.GRNEGI J,BA) 
GR(NEGI J)=ALF(FRNEGI J, GRNEG I J , F I NEG I J , G I NEG I J , BA > 

5 CONTINUE 
8887 K = K 
9999 KsK 

Ma6*ND+3 

C 

Ml = M-l 

DO 200 KM s 1 , Ml 

NST = 3 

M2 J = KM 

MKMsM-KM 

M2=MKM+1 

M21M s (M2-1)*M 

MM= ( M2-1 ) *M*M2 
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GRMM* GR < MM ) 

G I MM = GI(MM) 

[) a GRMM*GRMM + Si MM*G l MM 
01= l./D 

I F ( D*EPS ) 71 , 72, 72 
72 MONEY=0 

DO 2200 1=1, MKM 
!1M* ( I - 1 ) * M 
IM = I1M+M2 
IM*< I - 1 ) *M + M2 
GRIM = 6R< IM) 

G I I M = G I ( I M ) 

GAM 9 GRIM*GRMM * GIMM*3IIM 
GA= DI*GAM 

0E#DI*<GI IM*GRMM*GRIM*GIMM) 
FFICI)=GAM/D 

FFR(I) = (GIIM*GRMM » GRIM*GIMM)/D 
DE= FFR( I ) 

GA = FF I < I ) 

DO 200 J = 1, MKM 
D02200 J = 1# MKM 
IJ = I1M+J 
C ! Js ( I *1 ) *M+ J 

MJ= M21M**J 
C M J= ( M2-1 ) *M+J 

BE = GI(MJ) 

AL = GR(MJ) 

GR < 1 0 ) = GR(IJ)-GA*AU*D6*9E 
GI(IO) = G I ( I J ) * GA*8E-D£*At. 

2200 CONTINUE 
200 CONTINUE 
C REPLACE BETA(IU) 

C 

PI IN * G l ( 1 ) 

PIRN = QR(1) 

Ml = M-l 
G 

DO 103 K = 2» Ml 



KK»(K-1)*M+K 
GRKK ; GR(KK) 

G I KK= GI(KK) 

PJR = GRKK*PIRN-GIKK*PI IN' 

P I I s PI I N*GRKK + GIKK*PIRN 
PIRNsPIR 
103 PI IN»PI I 
C 

LJ=LJ-1 

RETURN 

C UNLIKELY EVENT A(MM)**2 ♦ B<MM)**2 TO SMALL 

71 NA=M-KM 
C 

DO 73 LL = 1 » NA 

LUVM=(LL-1)*M+J 
I F ( GR ( LUVM ) -DELTaX) 91,91,92 
91 I F C G I (LUVM)-DELTAX) 733,733,92 
73 CONTINUE 
C 

733 WRIT6(NP,1002)DELTAX,D,M2J 

1002 FORMAT ( 54H ****UNaBlE TO FIND AN ALPHA OR BETA LARGER THAN DELTA / 

1 21H IN SUBROUTINE DELY / 

2 10H DELTAXs , E20.8,17H OLD VALUE USEDs E20.8,10H COLUMN 
1 IS) 

I F ( D ) 72,999,72 

9? DO 77 J=1 , NA 
NA Js ( NA-1 ) *M + J 
K Js(LL-1)*M>J 

GR(NAJ) = GR(NAJ) * GR(KJ) 

G I (NAJ) = G I (NAJ) + G I < < J ) 

NANA = <NA-1)*M + NA 

D = GR(NANA)*GR(NANA)*GI (NANA)*GI (NANA) 

77 CONTINUE 
C 

4007 FORMAT (21H ADJUSTMENT ON COLUMN 1 5 > 

WRITE(NP,4007) M 
GO TO 72 

999 WRITE! NP, 1009) 

1009 FORMAT <1H050HNECESSARY TO ABORT DUE TQ SINGULARITY ) 

STOP 

DELTAX = 0,0 
L J=L J-l 

4 WRITE(NP,1005) 

RETURN 

END 
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SUBROUTINE UP ( PC# ND, NSTART, N* ) 

DIMENSION PC ( 18 ) 

C SUBROUTINE TO ESTABLISH CRITERION FOR CONVERGENCE 
C THE LAST THREE DETERMINANTS MUST BE MONOTONiC, 

EPSIL * ,075 
EP2 * 2,*EPSIL 
51 ND1 » ND-1 
ND2 * ND-2 

PA * PC(ND1)-PC(ND2) 

PD s PC(ND)-PC(ND1) 

P A65*PA 
PA76sPD 
P A 1 s ABS(PA) 

PDlsABS(PD) 

IFCPD1- PAU8-,55,55 
8 KsK 

IF(PA) 1,2,2 

1 PAi-1. 

GO TO 3 

2 P A = 1 . 

3 IF ( PD > 4,5,5 

4 PDs-1. 

GO TO 6 

5 PD=1, 

6 IF (PA+PD)7»55>7 

7 PDsPDl/PC(N01) 

PA*PA1/PC(ND2) 

PAsABS(PA) 

PDsABS(PD) 

IF (PA- EPSIL) 17,17,55 

17 IF (PD-EPS1L) 18,13,55 

18 KsK 

RAT J0=PA76/PA65 
IF (RATI0J55, 55,19 

19 IF (RATIO*. 8)54, 55, 55 

54 CALL PRE(ND,PC,NP> 

NSTART = 0 
RETURN 

55 CONTINUE 
NSTART s 1 
RETURN 
END 
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SUBROUTINE PRE ( ND , PC , NP ) 

DIMENSION PCI10) 

PREDICT THE CONVERGEED VALUES BASED ON THE LAST THREE 
DETERMINANTS 

REAL K2K1#K2K3,MU,KI,K2.K3 
DELTA*. 00001 
Cl = PC(ND-2) 

C2 = PC(ND-l) 

C3 8 PC(ND) 

KlsND-2 
K2 s ND-1 
K3 = ND 
K2K3* K2/K3 
K2K1*K2/K1 

WRITTEN FOR THE INFINITE DETERMINANT SUBROUTINE WITH ASYMPTOTIC 
LIMITS. 

PC ( 2 ) = 0 

MU = (C1-C2)/(C2-C3> 


P = 

FP 
A 

THIS 
AM s 
AP s 
M20*l 
M100*100 
Ml s 1 
13 K s K 


DO 

1 

M s M2Q, 

AP 

s 

K2K1*AP 

AM 

s 

K2K3*AM 

FP 

s 

MU-MU* AM 


*P ) -K2i<l**P* 
OR POSITIVE 

Ml 0 0 , Ml 
“ AP*1 


, 00 001 

8 MU*(1.-K2K3* 
a ABS ( FP ) /FP 

gets negative 
1 
1 


1 

ONE 
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8 s ABS(FP)/FP 
1 F ( A + 8 ) 1.10,1 
C ZERO IS CHANG E OF SI3N 
1 CONTINUE 
PC ( 2 ) = 3 
RETURN 
11 KsK 

10 M50 = M + 50 

18 PN = M 

19 FA = MU*(l,-K2K3**P'l>-K2Kl**=N + l 
DO IS L = 1,50 

PNK2K3 * K2K3**PN 
PNK2K1= K2K1**Pn 

PN1 = PN +(MU*(1,-PNK2K3)-PN<2K1 + 1 ) / ( PNK2K1*AL0G ( K2K1 ) +MU*PNK2K3* 
1AL0G<K2K3) ) 

PN1FP=(PN1-PN)/PN1 
PN = PN1 

8 = MU*(l.-K2K3**PNl)-K2Kl**=Nl +1 
PC ( 3 ) =PN1 

IF(ABS(PN1FP)-UELTA)16,16,15 

15 FA s B 
14 PC ( 2 ) = 2 

RETURN 

16 PC(1) =(C2*K2Kl**PNl-Cl)/( K2K1**PN1-1. > 

PC ( 2 ) = 0 

RETURN 

END 
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SUBROUTINE S5A < OUT » 4K » I H , K , MORDER , F„ PB , P A , PC , lERR) 

INTEGER R21 
INTEGER Rl 
INTEGER R,R2,R22 

DIMENSION B!GA< 13,13) * A(13»13) »B( 13,13) «C(13*13) 

DIMENSION BIGAP(13>13«6) ,BP(13*13,6) ,CP(13»13,6) , AP(13,13,6) 
DIMENSION F(13),AC<13),AK(13) 

DIMENSION ARU3), AS ( 13 ) , P ( 13 ) f Q( 13 ) , PC( 13 ) »PA ( 13 ) , PB( 13 ) 


THE POLYNOMIAL COEFFICIENTS OF THE HIGHER ORDER SYSTEMS. 

THE 9 THE AND 12TH ORDER SYSTEMS HAVE ASSOCIATED 
WITH THEM CHARACTERISTIC POLYNOMIALS OF DEGREE 9 AND 12 
RESPECTIVELY, THE COEFFICIENTS OF THESE POLYNOMIALS ARE NEEDED TO 
DEFINE THE SHAR ACTE R I ST I C EQUATIONS OF THE SYSTEMS8 AND 
ARE COMPUIEO AS FOLLOWS,,. 

THE COEFFICIENTS DESIREQ ARE K1,,,,K9 

NP = 6 

DO A 5 M = l, K 
M2= 2*M 
M21=M2-1 
PBM2 = PB { M2 ) 

C0SB2M = C0SCPBM2) 

PAM2 = P A ( M2 ) 

E A2M = EXPI-PAM2) 

F2M = F ( M2 ) 

E2A2M = EXP(-2,*PAM2) 

F2M1=F(M21) 

AR ( M ) = 2.*EA2M*(F2M1*SIN(P8M2>+F2M*C0SB2M) 

A S ( M ) a -2,*F2M*E2A2M 
P ( M ) = -2.*EA2M*C0S32M 
0 ( M ) = E2A2M 

WRITE(NP,2000)P(M),3(M),AR(M),AS(M) 

45 CONTINUE 
KlsK*l 

DO 10 M*K1# MORDER 
M2 a 2*M 
M21 = M2-1 
PCM2 s -PC ( M2 ) 

EC2M = EXP ( PCM2 ) 

PCM21 =-PC(M21) 

EC2M1 * EXP(PCM21) 

F2M = F ( M2 ) 
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F2M1 = F(M21) 

AR(M) = 2.*(F2M1*EC2M1*F2M*EC2M) 

CC9 = PCM21+PCM2 
C CC9 = -PCM21-PCM2 

ECC9 = EXP ( CC9 > 

AS (M ) = -2,*ECC9*(F2M1+F2M) 

P ( M ) = - (EC2M1+EC2H) 

Q < H ) = ECC9 
WRITE(NP, 20000 )M 
20000 FORMAT ( 5H M = ,15/) 

WRITE(NP,2000)P(M),3(M),AR(M),AS(M) 

20 0 0 FORMAT(1H0,/15X,4hP(M) , 15 X , 4 HQ ( M ) , 15 X „4HR ( M ) , l5x , 4HS ( M ) , / 

1 /3X.4E19.6) 

10 CONTINUE 

NOW WE DEFINE A SUCCESSION OF TRIANGULAR ARRAYS, 

A(I,2R), ,,, 1= 1 . 2 , • , . > 2R R= 1.2, 3, 4 

ACCORDING TO Al,2) = P(l) 

A ( 2 , 2 ) s 0(1) 

A(I,2R> = BI3A(I.R),*Q(R) + B(I,R)*P(R) * CU.R) 

WHERE BIGAU.R) = A(I-2,2R-25 
B ( I,R)= A( I-1.2R-2) 

C( I » R ) = A ( I.2R-2) 

WHILE IN COMPUTATION TAKE,,. 

A ( I , 2R ) s 0 FOR I NEGATIVE 
AND TAKE A ( 0 , 2 R ) s 1 

IN ORDER TO HAVE A(0,2R> WE SHALL ADD 1 TO ALL I INDICES TO AVOID ZERO 
N0RDER=MORDER*2 +1 
NL1= NORDER +1 
DO 25 I s 1 , NL1 
DO 25 R = 1,NL1 
25 A ( I , R ) = 0, 

ADJUSTED TO AVOID ZERO INDICES A. 1.2) AND A(2,2) 

A ( 2 » 2 ) = 0(1) 

A ( 1 , 2 ) = P(l) 

R=MORDER 
MR2 = 2 . *MQRDER 
MR21= MR2+1 
DO 5 I = 1.MR2 
K I = l 

KI1 = KI-1 
K 12= KI-2 
DO 5 R = l.MORDER 
R2 = 2*R 
R22 = 2*R-2*1 

ADJUSTED TO AVOID ZERO INDICES 
R22=R2-1 
I F ( K I 2 ) 1,2,3 
1 A ( K ] 2 , R22 ) = 0 
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GO TO 7 

2 A(KI2,R22> s 1. 

7 IF(KU) 4,6,3 

4 A ( K I 1 , R22 ) * 0 
GO TO 3 

6 A ( K 1 1 , R22 ) = 1, 

3 B I GA ( I , R ) = A ( K 1 2 , R22 ) 

B ( I , R ) s A(KU,R22> 

C( I,R) s A(KJ,R22) 

R21 = R2+1 

A ( K J » R21 ) * B I G A ( I ,R>*Q(R> + B<I,R)*P(R) * C(I»R> 

5 CONTINUE 

A ( 3 , 3 ) = 0,0 

A ( 2 , 3 ) r 0,0 

A ( 1 , 3 ) = 0 
DO 21 R = 4 , R2 , 2 
R1 = R+l 
DO 21 I = 1,R2 
A ( I , R > = A ( 1 , R 1 ) 

21 A ( I , R 1 ) s 0 

C WRITE (NR, 300 ) ( ( ( I , J, A( I, J) ) , 1=1,13) , J=l,13> 

300 FORMAT ( IX, 2HA < , I5,H, I 5 , 2H ) =E20 , 6 ) 

NEXT, DEFINE THE SUCCESSION OF MODIFIED TRAlNGULAR ARRAYS 
A P ( I , 2R , M ) 1= 1,2,,,.,2R R = l,2, . . . MORDER-1 M = 1,2, . . .MQRDER 

NOTE THE SUCCESSION HALTS AT R=MORDER-l 

THE ARRAYS OF AP ARE FORMED FOR A GIVEN M EXACTLY AS THE A(I,R) 

EXCEPT BEFORE ST ART I NS REPLACE THE VALUE OF P(M) BY P(4) AND □ ( M > BY Q ( 4 > 
N0RDER=l3 

DO 30 I s 1, NORDER 
DO 30 R = 1, NORDER 
DO 30 M = l,MORDER 
30 A P ( 1 , R , M > = 0,0 

DO 58 M s l,MORDER 
TMPPsP(M) 

TMPO?Q(M) 

0(M)sQ(M0RDER) 

P ( M ) = P ( MQRDER ) 

AP(2,2,M) = Q(l) 

AP ( 1 , 2 , M ) = P(l) 

ML1 = MORDER -1 
MR 2 = 2*M0RDER 

CEXCEPT BEFORE STARTING REPLACE P<M> BY P(MORDER) AND SAME FOR Q ( M ) 

MR21 = MR2*1 
DO 50 I = 1 , MR2 
K I = I 

KIlsKM 

KI2sKI-2 
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C ADJUSTED TO AVOID ZERO INDICES 
DO 50 R * 1 , HL1 
R2 s 2*R 
R22 = R2*l 
I F ( K 1 2 ) 51/52/53 

51 AP(KI2,R22,M) = 0 
GO TO 57 

52 AP(KI2,R22/M) s 1 

57 IF(KIl) 54,56,53 
54 AP ( K 1 1 , R22 # M ) s 0 

GO TO 53 

56 AP(KU,R22,M) s 1, 

53 R I G AP ( I , R / M ) s AP ( K I 2 / R22 , M ) 

BP ( I,R,M> * AP ( K 1 1 / R22 , M ) 

CP( I/R.H) * AP( K I , R22 # M ) 

R21 * R2 + 1 

AP ( I / R21 1 M ) s 8 I GAP ( I , R # M ) *Q ( R > * BP ( I . R / M ) *P < R ) ♦ CP<I,R/M> 

50 CONTINUE 
Q ( M ) = TMPQ 
P ( M ) = TMPP 

58 CONTINUE 

C NOW BACK SHIFT 

Rs MORDER-1 
R2 = 2*R 

DO 121 I * 1/NORDER 
DO 121 Msl/MORDER 
DO 121 R = 4/R2,2 
R1 = R+l 

A P ( I,R,M> = AP ( I,R1,M) 

121 AP C I/Rl/M) s 0 
DO 122 M*l,MORDER 

AP(1,3,H)=0 
A P ( 2 / 3 , M ) = 0 

122 CONTINUE 

400 F0RMAT(2X,3HAP( , I5/1H, I5/1H, I5/3H)= £20. 6) 

IT IS THEN POSSIBLE TO COMPUTE THE COEFFICIENTS C(S) S = l» 2, . t . 8 
ACCORDING TO ... 

CONST 1 = 0 

IF(M0RDER-4) 199,200/201 
201 IF(MORDER-6)199,202,199 
200 DO 100 M = l/MORDER 
100 CONSTls CONSTl+AH(M) 

AC ( 1 > = A(l. 8) + CONST1 
CONST1 = 0 
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DO 110 M = 1,4 

110 CONST1 = CONST1 * ARCM) * AP(1,6,M) * AS(M) 

AC ( 2 ) = A ( 2 » 8 ) + CONST 1 
DO 1200 M5*3,7 
CONST 1 = 0 
DO 120 M*l, 4 
MSI = MS-1 
MS2 s MS-2 

CONST 1 = CONSTl * AR{M)*AP{MS1,6 ,M)+AS<M)bAP(MS2»6,M) 
120 CONTINUE 

AC ( MS ) = A ( MS , 8 ) + CONST1 

1200 CONTINUE 

C WRITE(NP,500X('1,AC(I)),I=1.8) 

500 FORMA T ( lHO , 1X,3HAC(,I5, 3H)= ,E20.6> 

CONSTl = 0 
DO 130 M = 1,4 

CONSTl = CONSTl * A3< M > *AP ( 6, 6, M > 

130 CONTINUE 

A C ( 8 ) s A ( 8 , 8 ) CONST! 

AK1 s EXP ( -PC ( 9 ) ) 

AKO = 2 . *T <9)*AKl 
A K { 1 ) = AC ( 1 > -AKl+AKO 
DO 150 MS = 2,8 
MSI = MS-1 

AK(MS) = AC ( MS ) - AK1 * AC ( MSI ) ♦ AK0*A(MS1,8) 

150 CONTINUE 

A K ( 9 ) = AKO * A ( 8 , 8 ) - AK 1 * AC(8) 

55 0 F0RMATI3X, / /, 3X4HaK0b,E 2O .6,5X,4HAKl*,E20 .6/) 

650 FORMAT ( lH , IX, /.2x, 4HAK( ,I 5 . 3 H)= E20.6) 
WRlTE(NP»650X(I,A<(m,I=l,9) 

JERR = 0 
RETURN 

202 CONSTl = 0 

DO 205 M=l,MORDER 
CONSTl = CONSTl ♦ A R ( M ) 

205 CONTINUE 

AK ( 1 ) = A ( 1 , 1 2 ) + CONSTl 
CONSTl = 0 

DO 210 M = 1, MORDER 

CONSTl = CONSTl * AR(M)*A?(l,lQ,M)+AS(M) 


91 


210 CONTINUE 

AK<2) * A < 2 , 12 ) ♦ C0NST1 

DO 215 HS=3> 11 

MSI * MS-1 

MS2 = MS-2 

CONST1 = 0 

DO 212 M» l.MORDtR 

CONS T 1 s CONST 1 * A R ( M ) * A P ( MSI , 1 0 , M ) ♦ AS ( M ) * AP ( MS2 , 10 , M ) 
212 CONTINUE 

AK(MS) = A ( MS * 12 ) ♦ C0NST1 
215 CONTINUE 
CONST1 = 0 
DO 220 M=l*MORDER 

CONST 1 = CONST1 ♦ A 3 ( M ) * AP<10,10#M> 

220 CONTINUE 

AK ( 12 > = A ( 12 » 1 2 ) ♦ CONST 1 
WRITE(NP,650 ) ( ( I » AK( I ) ) , 1=1, 12) 

I ERR* 0 
RETURN 

199 IERR=MORDfcR 
RETURN 
FND 
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SUBROUTINE SIMQ(A,N,Y) 

DIMENSION A < 1 ) , Y ( t ) , ICHS(15)»SV(15) 

C SUBROUTINE FOR SOL, V I N3 SIMULTANEOUS EQUATIONS USING KROUTS METHOD 

DO 10 00 1 =1 , N 
II = ( I -1 ) *N+ I 
SVC I ) = A( I I ) 

C SVC I )=A( I , I ) 

IF (SVC I ) >5,46,5 
5 Y ( I ) = Y ( I >/SV( I > 

DO 10 00 J = 1 > N 
I J = ( I-1)*N+J 

AC I J> = AC I J)/SV( I ) 

1000 CONTINUE 
C10 00 A ( I , J ) = A ( I , J> /SV ( I > 

DO 101 K=1»N 
KK = (K-1)*N+K 
C AMX = ABS ( A ( K, K ) ) 

AMX= ABS C A ( KK ) ) 

IMX = K 

DO 15 I =K » N 
IKs (I -1>*N+K 
IFCABSCAC IK) )-AMX) 15,15,14 
IF ( ABS ( A ( I , K ) ) - A MX ) 15,15,14 
14 AMX = ABS ( A ( I , K ) ) 

14 AMX = ABSC A( IK) ) 

I MX* I 

15 CONTINUE 
I F ( AMX >27,46,27 

46 N = -> N 
RETURN 

27 IF ( I MX*K ) 8 , 9 , 8 
8 DO 22 J=1,N 
KJ = <K-1)*N+J 
C TEMP=ACK,J> 

TEMP = ACKvl) 

IMXJ - C IMX-1)*N+J 
A ( K J ) = AC IMXJ) 

C A (K, J)=AC IMX, J) 
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C 22 A( I MX , J ) sTEMP 
22 A ( I MX J )sJEMP 
ICHG(K)=IMX 
TEMPbY(K) 

Y(K)s Y ( I MX ) 

Y { I MX ) s TEMP 
GO TO 10 
9 I CHG ( K ) =K 

C 10 A(K#K)=l,/A<K,K) 

10 A ( KK ) = 1,/A(KK) 

DO 33 J=1,N 
IF (J-K) 6,33,6 

6 KJ s (K-1)*N+J 

A ( K J ) = A(KJ)*A(KK) 

C 6 A(K,J)=(A(K,J) )*(A(K* K) ) 

33 CONTINUE 

Y ( K ) = Y ( K ) * A ( KK ) 

C Y ( K ) = Y(K)*A(K,K> 

DO 44 1=1, N 
IK = ( I -1 ) *N + K 

7 DO 45 J=1 , N 

I F C I -K >17, 44, 17 

17 IF ( K" J ) 18 » 45 , 18 

18 IJ = ( I-1)*N+J 
KJ s (K-1)*N+J 

A ( I J ) = A( I J) »(A( I < ) * A ( K J ) ) 

C 18 A(I,J)=A<I.J)-(A<I,K)>*<A<k,J)) 
45 CONTINUE 

C Y(I> s Y ( I ) - A ( I , K ) * Y < K ) 

Y ( I ) = YU >-A( I K ) * Y ( K ) 

44 CONTINUE 

DO 99 1=1, N 
I F ( I-K)26,99,26 
26 IK = ( I -1 ) *N+K 

AUK) = -A( JK)*A(KK) 
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C 26 A(I,K)s-(A<l,K))*<A(K,K)> 
99 CONTINUE 
101 CONTINUE 

DO 70 K = 1,N 
L a N*l -K 
K I = I CHG ( L ) 

IF (L-KI) 68,70,60 

68 DO 69 1 = 1, N 

IL 9 (I-1)*N*L 
C TEMP = A ( I , L ) 

TEMP = A ( 1 L ) 

I K I =( I-1)*N+KI 
C A( I,l> = A( I,KI ) 

A ( I L > = A ( I K I ) 

A C I K I ) = TEMP 
C A ( I » K I ) = TEMP 

69 CONTINUE 

70 CONTINUE 

DO 1001 1=1, N 
DO 1001 J=1,N 
I J = ( I-1)*N+J 
C1001 A ( I, J)=A( I, J)/SV(J) 

A( I J) = AC I J)/SV< J) 

1001 CONTINUE 
RETURN 
END 
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SUBROUTINE S6(0UT,SS,IN,K,KH4T.NP) 

DIMENSION K<i2>,KHATC12).S(12),SS(12> 

0 SUBROUTINE FOR DETERMINING THE COMMON POLYNOMIAL COEFFICIENTS 
REAL KHA T * K 

REAL KHATl, <HAT2, <HAT3#KHAT4.KHAT5»KHAT6#KHAT7,KHAT8 

1, KHAT9,KHAT10,KHATli,KMAT12» 

1 K1,K2,K3,K4,K5,K6» <7 , K8 , K9 , <1 0 , K1 1 , MU1 » MUlSO , MUMU , 

2 L4,L3,L5*MU3,MU2,MJU,MJUU,K12»MU2SQ 
NPs6 

KHAT1= KHAT(l) 

KHAT2=KKAT(2) 

KHAT3 = KHAT < 3 ) 

KHAT4 s KHATM) 

KHAT5 = KHA r<5) 

KHAT6=KHAT(6) 

KHAT7 = KHAT ( 7 ) 

KHAT8=KHAT(8) 

KHAT9=KHAT(9) 

KHATIOsKHAT (10) 

KHATU= KHATC11) 

KHAT12=KHAT(12) 

K 1 s K ( 1 ) 

K2« K(2) 

K3 * K ( 3 ) 

K4 = K ( 4 ) 

K5 * K ( 5 ) 

K6 * K ( 6 ) 

K 7 = K ( 7 ) 

K8*K ( 8 ) 

K 9 = K ( 9 ) 

MU1 = K9 / KHAT12 

MU2 = (K8~KHAT11*MU1)/KMAT12 

MU2SQ=MU2*MU2 

MU3 * <K7-KHAT10*MU1-KHAT11*MU2)/KHAT12 
MUlSO * MU1 * MU1 
MUMU = MU2/MU1SQ 
MUUU = MU2SQ/MU1SQ - MU3/MU1 
HU = KHATl - K 1 
HI 0 = KHAT2 - K2 
L4 * HI 0 -Kl*Hll 
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A12 = K7*Hll + K 8 + K5/MU1 - K6*MUMU - KHAT8 

L 5 = HI 1 

L3 = KHAT3 - K3 -K1*H10 ♦ Hll*<Kl*Kl - K2) 

A13 = K7* ( H10-K1*H11 ) ♦ KS*Wll + «9 + K6/MU1 - KHAT9 

A21 = K5 + K2/MJ1 - K3*MUMU *K4*MUUU/MUl -KHAT5 

A22 = K5*HH + K6 ♦ K3/MU1 - K4*MUMU - KHAT6 

A 23 = K5*<H10-K1*H11>*K&*H11+K7+K4/MU1-KHaT7 
A31 = K3 + 1/MJ1 - K1*MUMU+<K2*MUUU)/MU1 - KHAT3 

A32 = K3*HH + K4 * K1/MU1 - «2*MUMU- KHAT 4 

A33 = K3*(H10-K1<H11) ♦ K4*HU + K5 «■ K2/MU1 - KHAT5 

R1 = KHAT10 -K7*L3 - K8 • L4 - «9*L5 
B2 = KHAT8-K8-K5*l3-K6*L4-<7*L5 
B3 = KHAT6 - K6 - K3*L,3 - <4*L4 - K5*l5 
FROM CRAMERS RULE. . , 

D1 = DELTAl/D, D2 = DELTA2/D. 03 = DE|_TA3/D 
A2233 = A22* A33 - A32 * A23 
A 1233 = A 12 * A 3 3 - A32*A13 
A1223 = A12*A23 - A 22 *A13 
0 c A11*A2233 - A21*A1233 * A31*A1223 
I F < D ) 15,5,15 

1002 F0RMAT(1H1,46H***DE\10MI \JATOR D FOR CRAMERS RULE IS ZERO EOJ 
5 WRITE(NP#1002) 

STOP 

15 CONTINUE 

OELTAl = B1*A2233 - B2*A1233 * B3 *Al223 
A2133 = A21-A33 - A31*A23 
A 1 133 = A11*A33 - A 31 * A 1 3 
A 1 123 = All* A23 - A 21 * A 1 3 

DELTA2 = -BW2133 ♦ B2 * A1133 - B3 *A1123 
A2132 = A21*A32 - A31*A22 
A1132 = All * A32 - A 3 1 * A 12 
A1122 = All* A22 - A21*A12 

DELTA3 = Bl *A2l32 - B2*A1132 + 83 *All22 

D1 = DELTAl/D 
D2 = DELTA2/D 
03 = DELTA3/D 

C THE COEFFICIENTS THEN ARE GIVEN BY ... 

SS(1)=K1-D3 

SS < 2 ) = K2-D2-D3*SS(1) 

SS < 3 ) * K3-D1-D2*SS(1)-03*SS(2) 

SS ( 4 ) =K4 - 01*SS<1> - D2 *SS ( 2 > -D3*SS ( 3 ) 

SS < 5 ) = K5 - D1*SS(2) - D2*SS(3)-D3 *SS<4) 

SS ( 6 ) = K6 -D1*SS(3)-D2*SS{ 4>-D3*SS(5) 

RETURN 

END 


97 


o a 


♦DECK, TEA 

SUBROUTINE TEA( NO, N3, ROOTR, ROOTj.Xl, ETA) 

DIMENSION ROOTR(l),ROOTI(2),XI(l),ETA(l) 

SUBROUTINE FOR I HE GENERATION Fo THE PSI AND ETA VECTORS FOR EACH 
SET OF POLYNOMIAL ROOTS 
PI = 3,1415926 
P I TWO = PI/2, 

P I 2 = 2 . *P I 
P I ONE = l./PI 
P I NV = 1 . /P I 2 
DO 15 1=1, NO 

SQS = ROOTRC I )*R00TR( I I+ROOTI ( I )*ROOTl ( I > 

XI(I) = -PINV*AL0Q(S0S) 

I F ( ROOTR ( 1 ) >14,50,14 

14 GAMMA = AT AN ( ROOT I ( I > /ROOTR ( I ) > 

ETA ( I ) = GAMMA 
IF ( ROOTR ( I) > 10,50,15 

50 I F ( ROOT I ( 1 ) ) 51,52,53 

51 ETA(I) =-PITWO 
GO TO 15 

52 K = K 

53 E T A ( I ) = PITWO 
GO TO 15 

10 IF (ROOTI ( .1 > ) 30,40,35 
40 ETA(I) = 3,1415926 
GO TO 15 

35 E T A ( I ) = ET A ( 1 ) ♦ PI 
GO TO 15 

30 E T A C I >=ETA( I > - P I 

15 FT A ( 1 > = -PlONE*ETA( I > 

RETURN 
END 

♦ END 
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SUBROUT J Nfc PAT(A>3,X»Z»X9,Xl2,X6»D,E»E9,£l2,E6*G»MU»NF,RR9,RI9» 

1 RR12,RI12,RR6,Rl6, ><09, VJ D 1 ? , Y9 » Y12 , P l ROLO, P I RNEW . P I I OLD. Pi I NEW, 

2 S»T'MD19 4 NDU2,YI9,YI12) 

DIMENSION PI 10LU(1>,PHNEW<1>, S(6,14), T < 11 , 14 ) , ND I 9 ( 1 ) 

SUBROUTINE FOR TH£ GENERATION OF THE SUMMARY TABLES FOR THE NASA FLUTTER 
DIMENSION YI9(1)»YI12(1) 

DIMENSION P1R0LD<1>,PIRNEH(1),A<8,14),B(U,14),X9<1),X12<1),X6{1) 
DIMENSION D(l).fc9(l).E12(l).E6(l),RR9(l),Rl9(l),RR12<l).RU2(l) 
DIMENSION RR6 ( 1 ) ,RI6(1)#Y9(1).E(1),X(1),Z(1)»ND9(1).NU12(1),U2<1) 

REAL MU 
GAMMA=G 
ND = 6 
NP = 6 

WRI T E ( NP, 5 0 ) GAMMA ,MJ,NF,C<X(J)»D(J>,Y9<J>.YI9<J)),J=1,9),(<Z<J>, 

1 E< J) » Y12( J) . Yll2( Jl) , J=1 .12) 

WRITE(NPjS4) 

DO 25 J = 1,8 
L = ND9(J) 

WR I TE ( ND, 51 ) (<Y9(J),A(J.l+1).A(J,L-2),A(J#L-1>.A(J,L).ND9(J>>. 

1 PIRNEW(J)) 

WRITE(ND,55l) ( < Y I 9 ( J) # S ( J , L + 1 ) » S ( J * L “ 2 ) •S(J*L-1)*S(J»L) , N D 9 ( J > ), 

1 PI I NEW ( U ) ) 

25 CONTINUE 

54 FORMAT (55X»24HNJNTH ORDER SYSTEM /) 

55 FORMAT <55X,24HTwElFTH ORDER SYSTEM /) 

WR I TE ( NP, 55 ) 

DO 20 M=l,ll 
L=ND12(M) 

JM8= M 

WRITE(ND,51)Y12(JM8),8(JM8,L+1),B(JM8,L-2),B(JM8«L-1). 

1 8< JM8.L) ,N012(M) , 3 JRO u D< JM3) 

WR ITE(ND»55l) Y 112 ( JM8) , T( JM8.L+1 ) ,T ( JM8.L-2) ,T ( JM8.L-1) « 

1T{ JM8,L) ,ND12(M) ,PII0LD<JM8> 

20 CONTINUE 

WRJTE(NP,53) ( ( X 9 ( J ) , E 9 ( J ) * RR 9 ( J ) » R I 9 ( J ) ) » J = 1 » 9 ) , < ( XI? ( J ) * E 1 2 < J > , 

1RR12( J) ,R L12( J) ) > J=l, 12) , ( ( X6 ( J ) ,E6(J) ,RR6( J) , R 1 6 ( J ) ) > J = l, 5) 
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50 FORMAT ( 1H1>38X»55HSTA3I|_I TY OF DYNAMIC SYSTEMS WITH PERIODIC P ARAM 

IMETtRS / 

2 /4 9X , 4 Q HROTOR WITH FlAP AND LEAD-LAG HINGES 

3/55X, 7HGAMMA = G12.4/55X.7H MU = Gl2 , 4 /55X » 7H NF = 15/ 

452X.40HSINGULARITIES AND EVALUATION POINTS / 

5 3 0 X » * H XI ETA YR 

2 YI* / 55 X *N I 

6NTH ORDER SYSTEM* / 9 ( 3 0 X , 4E2 0 . 8/ ) /55X , * T WELF TH ORDER SYSTEM* 
4/12(30X,4E20,8/)/6UX»18HDETERMlNANT VALUES 
5/15X ,1HY,18X,7H0. =S T , 14 X , 7 HQElT A 1 * 1 4X , 7HDEL T A 2 

6 14X,7HUELTA 3.10X.6HND MAX 10H N /) 

51 FORMAT (1X»*REAL*5E2D ,8* 1 1 2 » F 1 0 *5) 

551 FORMAT (IX, *IMAG*5E2Q,3,2X,I10»F10.5) 

52 FORMAT(1X,5E20,8»2X, 1 15. Ft 0.5) 

53 FORMAT ( *l*/55Xi 21HCHARACTERI 5T IC VALUES//55X * 24HN I NTH ORDER SYSTEM 

2 //19X,2HXI ,19X,3HETA*9X,18HREAL PART OF ROOT *3X, 

4 19HIMAG, PARI OF ROOT /, 9 ( 10 X , 4E20 . 8 / ) » 

3/.55X.24HIWELFTH ORDER SYSTEM 

4//19X.2HXIil9X,3HETA9X,18HREAL PART OF ROOT ,3X, 

519HIMAG. PART OF ROOT /l2 ( 1 0 X , 4E20 . 8/ ) , 

6/55X.24HS1XTH ORDER SYSTEM 

7//19X,2HX 1 , 19X, 3HETA9X, 18HREAL PART OF ROOT ,3X, 

819HIMAG. PART OF ROOT /6 ( 1 0 X , 4E20 , 8 / ) , /1H1 ) 

RETURN 

END 

♦ END 
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APPENDIX C 


Method for Extraction of the Common 
Polynomial Factor From the Two 
Higher-Degree Polynomials 

Given the two polynomials 

N ( z) = z 9 + kiZ 8 + k2Z 7 + ... + k 9 z + k 9 
T(z) = z 12 + kiZ 11 + k2Z 10 + ... + kjiz + kj2 

where it is known a priori that N and T both contain a sixth-degree 
factor S, 

S ( z) = z 6 + a x z 5 + c^z 4 + ... + 05Z + a 6 . 


it is required to determine aj, ozr in terms of the coeffi- 

cients k and k . 

n n 

Note, first, that the ratio of T to N must reduce to the ratio 
of a sixth-degree polynomial to a cubic, the common sextic factor 
cancelling off: 


T ( z) 
N(z) 


z 8 + X xz 5 + 1 2 z ** + ... + XgZ + lg 

Z 3 + WxZ 2 + U 2 Z + 113 


(C- 1 ) 


Now, a set of nine linear algebraic equations can be formed, the 
solution of which gives the values of the X's and y's. This set 

is obtained by multiplying Eq. (C- 1 ) through by N(z) [z 3 + v^z 2 
+ y 2 z + J13] and grouping the coefficients of like powers of z. 
There is considerable redundancy, with a total of fourteen equa- 
tions in the nine unknowns. A convenient set of nine are (with 
the associated power of z from which each was obtained indicated 
on the left) 


101 



z 1 **: 



ki 

+ VI 

z 13 : 


A 

k 2 

A 

+ ki y 1 

+ y 2 

z 12 : 

k 3 

+ k2 y 1 

A 

+ k^ y 2 

+ V 3 

z 5 : 

kio + 

k 9 y j + 

kgy 2 + 

k 7 y 3 

z 7 : 

kg + 

k 7 y 1 + 

kg y 2 + 

k 5 V 3 

z 9 : 

kg + 

k 5 yi + 

k^y 2 + 

kg y 3 


z° : k 12 y 3 
z 1 5 kiiH3+k 12 P2 
z 2 : k 10 P 3 + kiiy 2 + k 12 yi 


= k i + X i 

== k2 + + A 2 

— k3 + k£X 1 + k^X2 + A3 

= kgXi + k 8 x 2 + k 7 A 3 + k 6 X 4 + ksXs 
+ ^4X5 

= kg + kyAj + kgX2 ^5A3 + k 4 X 4 
+ ^3 A 5 + ^2 X* g 

= kg + kgXj + k 4 X 2 kgXg + ^2A4 
+ 1 A 5 + A g 

= kgXg 

= kg X 6 + kg X 5 
= k 7 X 5 + k g X 5 + kgX4 


If the first three and last three of these equations are properly 
combined and substituted in the fourth, fifth and sixth equations, 
then a set of three equations for the coefficients y 1 , p 2 and y3 

is obtained: 


3 


I a . . y 

• 13 

3 = 1 


4 -j 


i = 1 , 2,3 ; (C- 2 ) 
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where 


an 
a 1 2 

“1 3 
“2 1 
“22 
“23 
“3 1 
“32 
“3 3 

Bl 

e 2 

33 


k 7 + 


*IL 

5 i 




k 7 h 


7 n l l 


+ k« 


+ It - ft- ' 


^ 7 C^i x o ” k'lhii] + k 8 hii + kg 4 - kg 


■ *■ + tt - *• + -^[(ft) 2 - ft] - ^ 


k 5 hi! 4 k 6 + -^ 2 - - k 4 - k 6 


kst^io* kihii] -l* kghix + k7 4 - ky 


= k 3 + ^ - k, f^r + ff[(f^] 2 - ff] - k 3 


k 3 hii + k 4 + —J- - k 2 - k 4 


k3 [hi o - kxhij] + k 4 hx i + ks + - k 


k i o - k 7 £ 3 - kg£ 4 - kg^-5 


ks - k 8 - k5-£. 3 - k 6 £t+ - k 7 £ 5 


^6 ~ ^6 - ^ 3^3 - k 4 £ 4 - k 5^5 
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while 


h 1 1 - k 1 “ k l 


h i o 



- k 2 


l 3 = k 3 - k 3 - kj (k 2 - k 2 ) + (ki - kj) (k 2 2 - k 2 ) 
= k 2 - k 2 - k x (kj - k x ) 

•£5 = k i “ k i 


and 


Cl 


k i2 


C2 - 


(k 8 - ki 1 % 1 ) 


■12 


?3 ~ “ k 7 “ k lo5l - ^~ i “ L C 2 

k 1 2 _ k 10 J 


Eqs. (c-2) are readily solved for m, y 2 and 11 3 , using Cramer's 
rule. The coefficients of the sextic now follow immediately, since 


S (z) 


N(z) 

z 3 + y x z 2 + y 2 z + p 3 
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from which 


01 = ki - V! 

O 2 = k 2 - Ji 2 " v l a i 

03 = k 3 - y 3 - y 2 0 1 “ Pi a 2 

04 = k 4 - y 30 1 - UZ°2 ~ y l 0 3 

05 = ^5 “ y3°2 “ 0 2 0 3 “ yi°y 

0 6 ” ^6 “ y3°3 ~ U2 a 4 " y 1 0 5 
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